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ABSTRACT 

Cosmological hydrodynamical simulations are a valuable tool for understanding the 
growth of large scale structure and the observables connected with this. Yet, comparably little 
attention has been given to validation studies of the properties of shocks and of the result- 
ing thermal gas between different numerical methods - something of immediate importance 
as gravitational shocks are responsible for generating most of the entropy of the large scale 
structure in the Universe. Here, we present results for the statistics of thermal gas and the 
shock wave properties for a large volume simulated with three different cosmological nu- 
merical codes: the Eulerian total variations diminishing code TVD, the Eulerian piecewise 
parabolic method-based code ENZO, and the Lagrangian smoothed-particle hydrodynamics 
code GADGET. Starting from a shared set of initial conditions, we present convergence tests 
for a cosmological volume of side-length 100Mpc//i, studying in detail the morphological 
and statistical properties of the thermal gas as a function of mass and spatial resolution in all 
codes. By applying shock finding methods to each code, we measure the statistics of shock 
waves and the related cosmic ray acceleration efficiencies, within the sample of simulations 
and for the results of the different approaches. We discuss the regimes of uncertainties and 
disagreement among codes, with a particular focus on the results at the scale of galaxy clus- 
ters. Even if the bulk of thermal and shock properties are reasonably in agreement among 
the three codes, yet some significant differences exist (especially between Eulerian methods 
and smoothed particle hydrodynamics). In particular, we report: a) differences of huge factors 
(~ 10 — 100) in the values of average gas density, temperature, entropy, Mach number and 
shock thermal energy flux in the most rarefied regions of the simulations (p/pcr < 1) between 
grid and SPH methods; b) the hint of an entropy core inside clusters simulated in grid codes; 
c) significantly different phase diagrams of shocked cells in grid codes compared to SPH; d) 
sizable differences in the morphologies of accretion shocks between grid and SPH methods. 

Key words: galaxy: clusters, general - methods: numerical - intergalactic medium - large- 
scale structure of Universe 



1 INTRODUCTION 

Cosmological numerical simulations are a pow^erful 
tool to investigate the properties of the Universe at the 
largest scales. From galaxy formation to the precise 
measurement of cosmological parameters, from the 
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propagation of ultra high cosmic rays to the grow^th of 
the non-thermal energy components of the intra clus- 
ter medium (e.g. magnetic field, relativistic particles), 
cosmological simulations represent an effective com- 
plement to theoretical models and observations (e.g. 
Borgani et al. 2008; Borgani & Kravtsov 2009 and 
Norman 2010 for recent reviews). In order to model 
the evolution of cosmic structures in the most reliable 
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way, numerical methods must follow the non-linear 
dynamics of the gas and dark matter (DM) assembly 
across a very large dynamical range (e.g. from scales 
of - (10^ - 10^)Mpc to - (1 - 10)kpc), over the age 
of the Universe. 

To accomplish this task, a number of finite differ- 
ence methods have been developed in the past, which 
can be broadly divided into 2 classes (e.g. Dolag et al. 
2008 for a modem review). "Lagrangian" methods dis- 
cretize baryon gas by mass, using a finite number of 
particles, and the equation of fluid-dynamics are solved 
with the approach of smoothed particle hydrodynamics 
(SPH, see Price 2008 and Springel 2010 for recent re- 
views). Further details of the SPH method investigated 
in this project will be discussed in Sec l2.3[ 

Contrarily, "Eulerian" methods discretize space, 
by dividing the computational domain into regular cells 
(with fixed or variable size), and the gas-dynamics 
is evolved by solving cell-to-cell interactions (e.g. 
Le Veque 1990 for a review). A variety of numer- 
ical schemes can be applied for the reconstruction 
of the gas velocity, density, and pressure fields for 
a given number of neighbors (e.g. piecewise linear 
method, Colella & Glaz 1985; piecewise parabolic 
method, Colella & Woodward 1984), as well as for 
the time integration of the fluxes across the cells (e.g. 
ROE method, Powell et al. 1999; HLL/HLLE method, 
Harten et al. 1983; HLLC method, Li 2005). Further 
details of the grid methods em p loy ed in this project 
will be presented in Sections UA\ and f2?2\ 

Despite the enormous progresses made since their 
first applications (e.g. Peebles 1978; Efstathiou & East- 
wood 1981; Davis et al. 1985; Efstathiou 1985;), the 
mutual convergence of the results of cosmological nu- 
merical methods is still matter of debate and research. 
This is true also for the most simple physical model- 
ing of large scale structures, where no forces other than 
gravity and pressure are taken into account. 

A few comparison works in the literature (e.g. 
Kang et al. 1994, Frenk et al.l999, O'Shea et al. 2005, 
Heitmann et al. 2008) have provided evidences that 
most of the relevant quantities involved in large scale 
structure dynamics are generally reproduced with sim- 
ilar accuracy by most codes on the market. The general 
findings suggest that the simplest clustering properties 
of DM, and their dependencies on assumed cosmolog- 
ical and numerical parameters are fairly well under- 
stood (e.g. Heitmann et al. 2008). 

A less satisfactory agreement is generally found 
when the properties of gas in different methods are 
compared, even when simple non-radiative numerical 
setups are considered. In simulations of galaxy clus- 
ters, for instance, the entropy profile, the baryon frac- 
tion and the X-ray luminosities are affected by the 
larger uncertainties among codes reaching differences 
up to a factor of a few (e.g. Frenk et al. 1999; O'Shea 
et al. 2005; Voit et al. 2005; Kravtsov et al. 2005; Ettori 
et al. 2006; Vazza et al. 2010). 

More recent works aiming at comparing different 
numerical methods in more idealized test cases (e.g. 
shock tubes, blast waves, halo profile stability, ram 



pressure stripping of substructure) produced additional 
insights in the ways in which the numerical implemen- 
tations of different codes work (e.g. Agertz et al. 2007, 
Tasker et al. 2008, Mitchell et al. 2009; Springel 2010; 
Robertson et al. 2010; Hess & Springel 2010; Merlin 
et al. 2010). One of the reported key findings is that the 
effective numerical viscosity acting within each code 
has a sizable impact on the overall evolution of quanti- 
ties tightly linked to ram pressure stripping, turbulence 
and shocks in the gas medium. 

Cosmological simulations also proved to be impor- 
tant tools to study the acceleration and evolution of 
cosmic ray particles (CR) in the Universe, and their 
connection to the observed statistics of non-thermal 
emission from galaxy clusters (e.g. Dolag et al. 2008 
for a review). Several mechanisms related to clus- 
ter mergers and to the accretion of matter can act as 
sources of non-thermal components in the ICM. The 
most important mechanism during cluster formation is 
likely diffusive shock acceleration (DSA): the thermal 
particles in the high-energy tail of the Maxwellian dis- 
tribution function are able to experience multiple scat- 
terings across the shock surface which can be modeled 
as a diffusion process. This leads to an exponential gain 
of energy and an exponential loss of the number of par- 
ticles which results in a power-law distribution in par- 
ticle momentum extending into the relativistic regime 
and giving rise to so-called cosmic rays (e.g. Bell 1978; 
Blandford & Ostriker 1978; Drury & Volk 1981; see 
also Kang & Ryu 2010 and Caprioli et al. 2010 for re- 
cent reviews). 

Energetic shocks generated by mergers are be- 
lieved to accelerate supra-thermal electrons from the 
thermal pool and explain the origin of radio relics 
(Ensslin et al. 1998; Rottiger et al. 1999; Pfrommer 
et al. 2008; Pfrommer 2008; Hoeft & Bniggen 2007; 
Battaglia et al. 2009; Skillman et al. 2011), while 
high-energy electrons accelerated at these shocks can 
produce X-rays and gamma-rays via inverse Comp- 
ton scattering off CMB photons (e.g., Sarazin 1999; 
Loeb & Waxmann 2000; Blasi 2001; Miniati 2003; 
Pfrommer et al. 2008; Pfrommer 2008). Relativis- 
tic hadrons accelerated at shocks can be advected 
in galaxy clusters and efficiently accumulated there 
(Volk, Aharonian, & Breitschwerdt 1996; Berezinsky, 
Blasi, & Ptuskin 1997), possibly leading to a sizable 
non-thermal component which could be detected by 
gamma-ray observations (e.g., Pfrommer & Ensslin 
2004; Blasi, Gabici & Brunetti 2007; Pfrommer et al. 
2007; Pfrommer 2008; Pinzke & Pfrommer 201 1). The 
re-acceleration of relativistic electrons by MHD turbu- 
lence can be responsible for the episodic diffuse radio 
emission observed in the form of radio halos (e.g., Pet- 
rosian & Bykov 2008; Brunetti et al. 2008; Brunetti 
& Lazarian 2011); in addition secondary particles in- 
jected in the ICM via proton-proton collisions may 
also produce detectable synchrotron radiation (e.g., 
Blasi & Colafrancesco 1999; Dolag & Enssil 2000; 
Miniati et al. 2001; Pfrommer et al. 2008; EnBlin et 
al.2011). 

The occurrence of shock waves in large scale struc- 
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tures has been studied in detail with cosmological nu- 
merical simulations (e.g., Miniati et al. 2001; Ryu 
et al. 2003; Pfrommer et al. 2006; Pfrommer et al. 
2007; Kang et al. 2007; Hoeft et al. 2008; Skillman 
et al. 2008; Vazza, Brunetti & Gheller 2009; Molnar 
et al.2009) or indirectly by the action of shock waves 
on radio plasma bubbles, employing a novel method 
of combining radio observations and analytical insight 
that is supported by idealized hydrodynamic simu- 
lations (e.g., EnBlin et al. 2001; Pfrommer & Jones 
201 1). Most of these numerical works agree on the fact 
that the bulk of the energy in the Universe is dissipated 
at relatively weak shocks, M ^ 2 — 3 (where M is the 
Mach number), internal to clusters, while strong and 
larger shocks are found outside of large scale struc- 
tures, M ^ 10 — 100, at the boundary layers between 
the "collapsing" and the "expanding" universe. How- 
ever when the properties of CR injection by DSA are 
compared across the different simulations, differences 
up to 1 — 2 orders of magnitude in various quantities are 
found, including (but not limited to) quantities such as 
the ratio of energies of CR and thermal gas, the spec- 
tral energy distribution (e.g., Miniati et al. 2002; Ryu 
et al. 2003; Pfrommer et al. 2006,2007,2008; Pfrom- 
mer 2008; Kang et al. 2007; Hoeft et al. 2008; Vazza 
et al. 2009,2010; Skillman et al. 2008,2011; Pinzke & 
Pfrommer 2011). This limits our present understand- 
ing of the main mechanism for the enrichment of CRs 
in the intra cluster medium. 

In this work, we explicitly aim at comparing three 
independent numerical approaches for cosmological 
simulations, applied to the evolution of a large volume 
of the Universe: the smoothed particle hydrodynam- 
ics code GADGET (Springel, Yoshida & White 2001; 
Springel 2005); the total variation diminishing code de- 
veloped by D. Ryu and collaborators (Ryu et al. 1993; 
Ryu et al. 2003) and the parabolic piecewise method 
ENZO, developed by G.Bryan and collaborators (e.g. 
O'Shea et al. 2004; Norman et al. 2007). 

We adopted a set of shared identical initial con- 
ditions generated at different resolution, and we re- 
simulated them with the three codes; the output of all 
runs were then compared in detail, looking at the con- 
vergence of several thermal and non-thermal properties 
across the various codes and for different numerical 
resolutions. We chose the simplest possible physical 
setup for this project, and include only non-radiative 
physics (i.e., no radiative cooling, no UV radiation 
background from primordial stars, no magnetic fields, 
etc.). 

This approach helps us to understand which dif- 
ferences are due to the numerical methods (e.g. "La- 
grangian" versus "Eulerian" method for gas dynam- 
ics) and which are due to the post-processing (e.g. 
temperature-based method to detect shocks versus ve- 
locity based methods). Also, this approach helps in as- 
sessing some of the more robust findings of present 
cosmological simulations, and determines the mini- 
mum resolution requirements needed to achieve a good 
convergence independent of the particular adopted nu- 
merical method. 



Table 1. Details of the simulations run for this comparison project. First 
column: name of the run; second column: initial redshift of the simula- 
tion; third column: mass resolution for Dark Matter particles; fourth col- 
umn: softening length (for SPH runs) or uniform mesh spacing (for ENZO 
and TVD) employed in the runs. 



GADGET 



Run 




Mdm [Mo/Zz] 


Rsoft [kpc//i] 


o4 


34. o3 


2.4 • 10 


31. U 


128 


44.77 


3.0 • lO^o 


15.75 


256 


55.92 


3.76 • 10^ 


7.875 


ENZO 








Run 


^in 




Ax [kpc//i] 


64 


34.63 


2.4- lO^i 


1562.5 


128 


44.77 


3.0- lO^o 


781.25 


256 


55.92 


3.76 • 10^ 


390.625 


512 


67.99 


4.7- 10^ 


195.31 


TVD 








Run 




Mdm [Mq/K] 


Ax [kpc//i] 


64-32 


34.63 


3..0- 10^2 


1562.5 


128-64 


44.77 


2.4- lO^i 


781.25 


256-128 


55.92 


3.0 • lO^o 


390.625 


512-256 


67.99 


3.76 • 10^ 


195.31 



The paper is organized as follows. In Sec|2l we 
give a brief description on the underlying numerical 
schemes of these codes. Based on the simulations us- 
ing different resolutions we present a comparison of 
the general distribution statistics for dark matter and 
thermal gas in sections |4] and O In particular, we fo- 
cus on the galaxy clust ers p roperties according to the 
various codes in section [531 and present an exploratory 
test showing important differences between the under- 
lying numerical, hydrodynamical schemes (specially 
between PPM and SPH) in th e matter accretion pat- 
tern inside halos in section [54l We then apply different 
shock detecting schemes in section 16.11 to the various 
re-simulations and we present results for the ch arac - 
teriz ation of shock waves in all codes in sections (63]- 
16.51 We particularly focus on shocks in galaxy clusters, 
their properties and their role in the acceleration of CR 
predicted according to the d ifferent, underlying numer- 
ical schemes in section [6^ 



2 NUMERICAL CODES 

2.1 Eulerian method: ENZO PPM 

ENZO is an adaptive mesh refinement (AMR) cosmo- 
logical hybrid code originally written by Greg Bryan 
and Michael Norman (Bryan & Norman 1997, 1998; 
O'Shea et al. 2004; Norman et al. 2007). It couples a 
particle-mesh solver with an adaptive mesh method for 
ideal fluid-dynamics (Berger & Colella, 1989). 

ENZO uses a particle-mesh N-body method (PM) 
to follow the dynamics of collision-less systems. This 
method computes trajectories of a representative sam- 
ple of individual DM particles and it is much more ef- 
ficient than a direct solution of the Boltzmann equation 
in most astrophysical situations. 

DM particles are distributed onto a regular grid 
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using the cloud-in-cell (CIC) interpolation technique, 
forming a spatially discretized DM density field. After 
sampling dark matter density onto the grid and adding 
baryon density (calculated in the hydro method of the 
code), the gravitational potential is calculated on the 
periodic root grid using Fast Fourier Transform algo- 
rithms, and finally solving the elliptic Poisson's equa- 
tion. 

The effective force resolution of a PM calculation 
is approximately twice as coarse as the grid spacing at a 
given level of resolution. The potential is solved in each 
grid cell; however, the quantity of interest, namely the 
acceleration, is the gradient of the potential, and hence 
two potential values are required to calculate this. 

In the case of ENZO simulations employing AMR, 
the potential is recursively computed within sub-grids 
at a higher resolution and the boundary conditions are 
interpolated from the potential values of the parent 
grid. Then a multi-grid relaxation technique is adopted 
to compute the gravitational force for each cell within 
sub-grids (e.g. O'Shea et al. 2006). This enables the use 
of a gravitational softening of the order of the highest 
resolution available in the simulation; however in this 
project we did not use AMR capabilities of ENZO, and 
thus the maximum available softening is the fixed res- 
olution of the adopted mesh. 

As hydrodynamical solver, ENZO adopts the Eule- 
rian Piecewise Parabolic Method (PPM, Woodward & 
Colella, 1984). The PPM algorithm belongs to a class 
of schemes in which an accurate representation of flow 
discontinuities is made possible by building into the 
numerical method the calculation of the propagation 
and interaction of non-linear waves. It is a higher or- 
der extension of Godunov's shock capturing method 
(Godunov 1959). It is at least second-order accurate 
in space (up to the fourth-order, in the case of smooth 
flows and small time-steps) and second-order accurate 
in time. This leads to an optimal treatment of energy 
conversion processes, to the minimization of errors due 
to the finite size of the cells of the grid and to a spatial 
resolution close to the nominal one. 

In order to treat more accurately bulk hypersonic 
motions, where the kinetic energy of the gas can dom- 
inate the internal energy by many orders of magnitude, 
both the gas internal energy equation and total energy 
equation are solved everywhere on the grid at all times. 
This dual energy formulation ensures that the method 
produces the correct entropy jump at strong shocks and 
also yields accurate pressures and temperatures in cos- 
mological hypersonic flows. 

This works uses the public 1.0.1 version of 
ENZOQ. To simplify the comparison with the other 
codes of this project, this work employs a fixed grid 
only instead of the adaptive multilevel grids and ad- 
ditional physics (e.g. star formation, re-ionization and 
cooling processes) which are powerful tools in ENZO. 



http://lca.ucsd.edu/software/enzo/vl .0. 1/download/ 



2.2 Eulerian method: Cosmological TVD 

The cosmological code created by Ryu et al. (1993) 
is based on the Harten (1983) Total Variation Dimin- 
ishing (TVD) scheme. It is a flux-based Eulerian code 
with second-order accuracy in space and time. It cap- 
tures shocks within two to three cells without gener- 
ating oscillations, but limiting the numerical flux ac- 
cording to the TVD scheme instead of adding a sim- 
ple artificial viscosity. Several important improvements 
were made while incorporating the TVD scheme into 
the cosmological code. The numerical artificial heat- 
ing around the extremely supersonic flows where the 
bulk kinetic energy is much greater than the thermal 
energy is reduced; this was achieved by following the 
adiabatic changes of the thermal energy using a mod- 
ified -entropy equation instead of using the total en- 
ergy equation. The leakage of the gravitational energy 
into the thermal energy in regions of supersonic flows 
was prevented by including the effects of the gravita- 
tional force only to the momentum and kinetic energy 
and keeping the thermal energy rather than solving the 
conservation of the total energy. Also, a correction due 
to the mass diffusion under the gravitational field has 
been added in the gravitational force term in order to 
obtain better conservation of the total energy and to 
satisfy the cosmic energy equation. Additional details 
can be found in Ryu et al. (1993) and Ryu et al. (2003). 

The treatment of gravity and DM particle dynam- 
ics follows the Particl e M esh approach on a fixed res- 
olution grid (see Sec l2.1l) . Additionally, in this code 
there is the possibility of using a number of DM parti- 
cles smaller than the total number of cells in the grid, 
in order to spare memory usag e. T his is motivated by 
the fact that, as stressed in Sec l2.1l in the PM scheme 
the effective force resolution is approximately twice as 
coarse as the mesh spacing. Therefore, adopting a num- 
ber of DM particles which is for a grid, has 
a very little or negligible difference in the final accu- 
racy of the resulting potential and accelerations. 



2.3 Smoothed Particle Hydrodynamics: GADGETS 

We compare Eulerian methods with the parallel 
TreeSPH code GADGET3 (Springel 2005), which 
combines smoothed particle hydrodynamics with a hi- 
erarchical TreePM algorithm for gravitational forces. 
SPH uses a set of tracer particles to discretize mass el- 
ements of the fluid. Continuous fluid quantities are es- 
timated by a kernel interpolation technique (e.g. Mon- 
aghan 1992). The equation of motion for these tracer 
particles can be derived (by applying the variational 
principle) from the Lagrangian of such system. The 
thermodynamic state of each fluid element may ei- 
ther be defined in terms of its thermal energy per unit 
mass, Ui, or in terms of the entropy per unit mass, Si. 
The latter is used as the independent thermodynamic 
variable evolved in SPH, as discussed in full detail by 
Springel & Hemquist (2002). The adaptive smoothing 
lengths hi of each SPH particle are defined such that 
their kernel volumes contain a constant mass for the 
estimated density (e.g. corresponding to the mass of 
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TV = 64 particles is a common choices). Accounting 
for the fact that then the adaptive smoothing lengths 
hi are a function of density allows SPH to be formu- 
lated so that both energy and entropy are manifestly 
conserved. Provided there are no shocks and no ex- 
ternal sources of heat, the derivation of equations for 
the reversible fluid dynamics in SPH is straightforward 
(see Price 2008 and Springel 2010 for recent reviews 
on SPH). However, flows of ideal gases can easily de- 
velop discontinuities where entropy must be generated 
by micro-physics. Such shocks need to be captured by 
an artificial viscosity in SPH, which is active only when 
fluid elements approach one another in space, prevent- 
ing particle interpenetration and transforming kinetic 
energy irreversibly into heat (e.g. Monaghan & Gin- 
gold 1983). Modem schemes like GADGETS make 
also use of an artificial viscosity based on an analogy 
with Riemann solutions of compressible gas dynamics, 
as proposed by Monaghan 1997; additional viscosity- 
limiters are also introduced in GADGETS in the pres- 
ence of strong shear flows to alleviate spurious angular 
momentum transport (Steinmetz 1996). 

Both the collision-less dark matter and the gaseous 
fluid are represented by particles, allowing the self- 
gravity of both components to be computed with grav- 
itational N-body methods. GADGETS allows the pure 
tree algorithm to be replaced by a hybrid method con- 
sisting of a synthesis of the particle-mesh method and 
the tree algorithm, with significant reduction of the 
computational effort. 

The effective force resolution is controlled by the 
gravitational softening i?soft used in the tree part as 
listed in the last column of table 1 for the different sim- 
ulations and the particles are allowed to have individ- 
ual time steps, based on different time stepping criteria 
(see Springel 2005 for details). 



3 INITIAL CONDITIONS 

We have assumed a "concordance" model, with density 
parameters ricurv = 0, = 0.043, ^dm 0.227, 
= 0.73, Hubble parameter h = 0.70, a power 
spectrum with slope n = 1 and a normalization of 
the primordial matter power spectrum erg = 1.2. The 
(78 parameter is intentionally set to a larger value com- 
pared to recent estimate from CMB data (e.g. Spergel 
et al. 2007) in order to enhance the probability of form- 
ing massive halos within the simulated volume of side 
100Mpc//i. Any modeling of cooling, radiative and 
heating processes for the gas component is neglected, 
and therefore the thermal history of cosmic gas here is 
mainly driven by shock waves induced by gravity. Ta- 
ble [T] lists the main parameters of all simulations run 
for the project. 

The initial displacements and velocities of DM par- 
ticles were identical for all codes; the numbers of DM 
particles adopted are 512^, 256^, 128^ and 64^. The 
GADGETS simulations preliminary looked remark- 
ably converged with resolution already at 256^, and 
therefore we choose to skip the production of the 512^ 



case in SPH, in order to spare computational resources 
time. The initial redshift of simulations were computed 
in order to reach the same growth rate at z = for the 
smallest available density perturbations: zin = 67.99, 
Zin = 55.92, Zin = 44.77 and zin = 34.63 for the 
different resolutions, respectivel}0. 

Usually in SPH cosmological runs both the DM 
and the gas particle distributions are perturbed in 
their initial positions and velocities according to the 
Zel'Dovich approximation (e.g. Dolag et al. 2008 for a 
review). In grid runs, on the other hand, the initial gas 
distribution is at rest compared to the DM initial veloc- 
ities. Since computing exactly the same initial pertur- 
bation in velocity for SPH particles and cells is not a 
trivial issue, for the sake of simplicity in this work we 
neglected initial perturbations in velocities also for the 
SPH distribution. 

In the following we will refer to a given run ac- 
cording to the number of its gas particles or gas cells; 
in the case of the TVD code, the number of DM par- 
ticles is kept 8 tim es smaller than the number of gas 
cells (see Sec i2.2l) . In what follows, we will typically 
refer to "^e//^convergence" meaning the convergence 
of a code with respect to increasing resolution, and 
to "cr6>5'5'-convergence" meaning the convergence be- 
tween different codes, at a given resolution. 



4 DARK MATTER PROPERTIES 

A number of works in the literature have shown that 
present day numerical codes at their best achieve an 
agreement within ^5 — 10 per cent on the mass func- 
tions of halos (e.g. Frenk et al. 1999; O'Shea et al. 
2005; Heitmann et al. 2007). However, subtle differ- 
ences in the adopted numerical methods should be re- 
sponsible for the exact shape of the inner DM profiles 
(e.g. Bullock et al. 2001; Warren et al. 2006). 

We compared the properties of the DM component 
for all resolutions and codes in order to ensure that the 
distribution of DM in our simulated large scale struc- 
tures is characterized by a similar degree of intrinsic 
"scatter" reported in the literature. 

The most important statistics related to DM is the 
mass function of halos, for which analytical solutions 
as a function of cosmological parameters are avail- 
able (e.g. Press & Schecter 1974; Sheth & Tormen 
1999). We report in Fig. [T] the cumulative mass func- 
tions (DM plus gas) for all runs in the project. The 
virial mass, Myir, is customarily defined as the spher- 
ical over-density of gas+DM, enclosing a mean over- 
density of ^ 109pcr, where pcr ^ 9.31 • 10~^^^/cm^ 
is the critical density of the universe (e.g. Eke et al. 
1998). The virial halo masses are computed using the 
same halo finder in all codes, based on the gas+DM 
spherical over-density. In the case of grid runs, the cells 
distributions have been converted into a distribution of 



^ The initial conditions used in this Project are pubHc and accessible at: 
|http ://canopus .cnu. ac .kr/shocks/caseO/ 1 
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Figure 1. Cumulative mass functions of the virialised halos in the various runs. In both panels the GADGET results are reported in black for the various 
resolutions, while the left panel reports the mass functions from ENZO runs (in blue) and the right panel reports the mass functions from TVD runs (in red). 
The Sheth & Tormen (1999) mass function is shown for reference in bold (grey lines), with the thin lines showing the Poisson errors. The vertical lines indicate 
the minimum mass resolution for each cluster run, as outlined in Sec|4] 




Figure 2. Baryon fraction for all halos in the three codes, at all resolutions. The vertical lines mark the minimum mass resolution criterion outlined in SeclH 




B[H^] 0,01 0.10 1.00 10.00 



Figure 3. Radial profile of DM density for the most massive galaxy cluster in our volume, for the three codes; the profiles of the 256 run in GADGETS are 
reported for comparison in the last two panels. The vertical lines in each panels show the value of gravitational softening for each run. 
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Figure 4. Maps of projected mass-weighted temperature through the full cosmological volume of 100Mpc//i, for the three most resolved runs of our project. 




Figure 5. Maps of projected mass-weighted temperature for a sub region with the side of 40Mpc//i for all codes and resolutions. 
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particles, in order to apply exactly the same procedure 
used to analyse GADGET runs. 

In order to compare different codes and resolu- 
tion, it is useful to assign a "formal" resolution to each 
run. This allows us to understand which halos in our 
simulations are suitable for "convergence" studies and 
which are not, because of under-sampling problems at 
a given cluster size. Even if in GADGET runs the mass 
functions are resolved down to the smallest halos (with 
< 20 particles). Power et al. (2003) showed that con- 
vergence in the inner dynamical structures of halos is 
achieved with at least ^ 500 particles inside i?vii0- 

We preliminary consider that the resolution limit 
in GADGET is achieved with 500 DM particles within 
the virial radius. For grid runs, we apply the follow- 
ing empirical approach: we consider only halos whose 
virial radius is resolved with at least 500 cells, and we 
assign ^formal minimum mass to have halos "suitable 
for convergence" taking the corresponding virial mass, 
extracted from the theoretical i?vir versus Myir rela- 
tion. 

The corresponding minimum masses for all codes 
and resolutions are shown as vertical lines in Fig. [TJ 
Although this methods is rather artificial, we find it 
predicts rather well the convergence observed for ha- 
los in grid codes, which takes place at larger masses 
compared to corresponding GADGET runs at the same 
DM mass resolution. For instance, GADGET run 256 
shows a halos mass function which is converged down 
to masses of 2 • lO^^M©, while run 256 in ENZO and 
TVD achieve convergence only for halos with masses 
larger than 5 - lO^^M©. 

Therefore, we would expect to see cross- 
convergence of the virial parameter for none but the 
largest halos in grid results, while we expect good self- 
convergence across a larger range of masses in GAD- 
GET. 

A similar trend is also observed in the baryon frac- 
tion of halos in the various run, as reported in Fig. [2l 
The baryon fraction in GADGET is rather perfectly 
converged at all resolutions for M > 10^^ M^/h, with 
a value of ft ^ 0.9/cos, where /cos = ^h/{^h + 
f^dm) = 0.159 is the cosmic baryon fraction in our 
runs) . In grid codes, the convergence to a slightly 
larger baryon fraction, /b ^ 0.95/cos. seems to be 
reached only for masses larger than M > IO^^Mq/Zi 
(as for the halos mass functions, ENZO shows the a 
slower rate of convergence compared to TVD). 

The radial profiles of DM mass density for the 
most massive galaxy cluster in our sample are shown 
in Fig. [3] for various resolutions. All profiles in GAD- 
GETS runs are remarkably self converged, while the 
profiles of DM in both grid methods present a slower 
rate of convergence. At the best available resolution, 
the grid codes agree at the percent level with the ref- 
erence profile of GADGETS runs, with sizable differ- 

^ We notice, however, that tests with radiative runs have shown that a larger 
number of particles, ~ 1000 — 5000, may required to achieve a good 
convergence in the X-ray luminosities of clusters (e.g. Valdarnini, Ghizzardi 
& Bonometto 1999; Valdarnini 2002). 




Figure 6. One dimensional distribution of gas mass density (lower lines) 
and volume- weighted temperature (upper lines) for a line crossing our sim- 
ulated volume, for all 256^ runs. 



ences only in the core region of the cluster, < 0.1i?vir, 
due to the well known lack of force resolution in the 
PM method (the softening length for the gravity force 
in the 512^ runs is 293kpc). 

Overall, the trend found are in line with those re- 
ported by O'Shea et al. (2005) and Heitmann et al. 
(2008). Based on our results, we suggest that the repre- 
sentation of the underlying DM distribution is similar 
to what can be found in the recent literature, and that 
the bulk of differences that will be reported in the next 
Sections are mostly connected with a different model- 
ing of hydrodynamics in the various methods. 



5 BARYONIC MATTER PROPERTIES 

In the following Sections we compare the distributions 
of several gas thermodynamical variables in all runs 
as a function of numerical resolution. The final goal 
is to identify which are the cosmic environments and 
minimum resolution requirements necessary to achieve 
a good convergence in the estimates provided by the 
different methods. 



5.1 Maps 

A preliminary inspection of the morphological distri- 
bution of baryon gas in the cosmic structures captured 
by all methods ensure that at a zero order all simu- 
lations correctly sample a cosmological volume with 
identical density fluctuations. In Figure[6]we report the 
one-dimensional behavior of gas density and gas tem- 
perature along a line crossing the position of the most 
massive galaxy cluster in the volume, for all 256^ runs. 
The spatial distribution of gas density is well matched 
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in all codes, and in particular the positions of the gas 
density peaks associated with halos and filaments agree 
within a 1 — 2 cells accuracy (i.e. ^ 400 — 800kpc//i at 
this resolution). The one-dimensional gas temperature 
profiles show very similar maxima near the gas density 
peaks, but sizable differences can be found in the outer 
regions. The bulk of the difference here is however a 
simple effect of the variable smoothing length in GAD- 
GETS, which provide a coarser resolution compared to 
grid codes for the regions outside of clusters. 

In the panels in Fig. |4] we report the maps of 
projected mass-weighted temperature across the simu- 
lated volume, for the most resolved runs of the project 
(run 256 for GADGETS and runs 512 for ENZO and 
TVD). The trend with resolution of the projected mass- 
weighted temperature, at all resolutions, is reported in 
Fig.oJfor a sub-volume of 40Mpc//i inside the cosmo- 
logical box. 

To readily compare Lagrangian and Eulerian data 
at the same spatial resolution, the gas fields of GAD- 
GETS runs have been interpolated onto a regular grid, 
with resolution equal to that of the corresponding grid 
runs, using the same SPH kernel employed during the 
simulation for each gas particle. 

In GADGETS, over dense non-linear structures 
(e.g. halos and sub-halos) are very similarly recon- 
structed at all resolutions, while structures at about the 
critical density (e.g. cosmic filament) start being re- 
solved only at sufficiently high DM mass resolution. 
The opposite trend appear in grid codes, where large 
scale patterns are soon reconstructed at all resolutions, 
while a clear modeling of the smaller halos and clus- 
ter satellites is achieved only approaching the highest 
available resolutions. 



5.2 Distribution Functions 

A quantitative analysis of the differences between the 
codes is performed by studying the volume-weighted 
distribution functions of gas density and gas tempera- 
ture at increasing resolution, as shown in Fig. Ul Fig- 
ure [8] further shows the cross-comparison between the 
highest resolution runs available for each code. 

In this case, we adopt volume-weighted statistics 
for each bin in gas density/temperature. Despite the 
obvious fact that volume-weighted distributions can- 
not be translated into observable quantities (since the 
convolution of the two does not provide the total gas 
energy within the simulated volume) we find this ap- 
proach useful to focus on the properties of the low 
density, volume-filling baryon gas around large scale 
structures. Our purpose here is to highlight the differ- 
ences in the modeling of the lower density baryon gas 
at large scales (which encompasses filaments and clus- 
ters of galaxies) in the different numerical methods. 
This can also be readily compared with the early com- 
parison work of Kang et al. (1994). In addition, these 
volume filling regions are expected to be an important 
site of acceleration of relativistic particles, via direct 
shock acceleration at strong shocks (e.g. Miniati et al. 
2001; Ryu et al. 200S; Pfrommer et al. 2006; Vazza, 
Brunetti & Gheller 2009). 
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Figure 8. Cross convergence of volume weighted gas density and gas tem- 
perature distributions for GADGETS run with 256^ and grid runs with 
5123. 



In the following Section (Sect. 15. SI) we will rather 
refer to mass-weighted profiles of gas density and gas 
temperature, since they are closely related to the ther- 
malization properties of internal merger shocks inside 
clusters. 

As expected, the cross-convergence between dif- 
ferent codes is more satisfactory when resolution is 
increased: the density distributions runs with > 256^ 
DM particles (i.e. with rridm ^ 4.5 • 10^ M^/h), have 
the same average value in all codes within a 20 — 30 
per cent scatter. The largest and the smallest gas densi- 
ties are similar within a factor of ^ 2, and GADGETS 
produces the most extreme values in both cases. GAD- 
GETS runs are also the ones which provide the largest 
degree of self-convergence, with very similar outputs 
at all investigated resolutions. 

In the case of temperature distributions, ENZO 
presents the larger degree of self-convergence (within 
a factor of ^ 10 per cent) at all resolutions, while the 
other codes show significant evolution with resolution, 
especially at temperatures below T < 10^~^ K. 

We note that different floors in the value of tem- 
perature were adopted in the three codes, to limit the 
lowest temperature available to cells/particles of the 
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Figure 7. Gas density and volume-weighted gas temperature distributions for all resolutions and all codes. The critical density of baryons is pcr,b ~ 4.0 ■ 



simulated volume. For each code we used the temper- 
ature floor usually adopted by each simulator: a mini- 
mum temperature of Tq = 1 K in allowed in ENZO, 
To = 2 K in TVD and To = 24 K in GADGETS. 
This explains the different piling of cells/particles in 
the temperature distributions below T < 50 we 
also made sure that the adoption of different floor in 
temperature does not affect in any way the temperature 
distribution above the adopted Tc0 . 

On the other hand, the temperature distributions 
found in the simulations become quite similar for T > 
10^ K, which would correspond to the typical virial 
temperatures of collapsed halos; this is in line with the 
early findings reported by Kang et al. (1994), and later 
by O' Shea et al. (2005). 

We conclude that even if the gas mass distribution 
within halos is rather convergent in all codes for (for 
a DM mass resolution of rridm ^ 4.5 • lO^MQ/h), 
the convergence in the gas temperature distribution is 
generally not yet reached, and the cross-convergence 
between codes is not achieved for all regions where 



^ It should be stressed that all most recent simulations model the action 
of the re-ionization background, hence increasing the minimum tempera- 
ture in the simulations to much larger values, To ~ 10^ — 10^ K (e.g. 
Vazza, Brunetti & Gheller 2009). Therefore the analysis of the temperature 
distribution we present here is meant to pinpoint the numerical problems 
of the various methods, while the differences between runs employing re- 
ionization would be much smaller. 



T < 10^ — 10^, for the resolutions investigated in this 
project. 

In these regimes, some amount of spurious numer- 
ical heating can be expected due to the graininess of 
DM mass distributions, which makes two-body heat- 
ing a likely channel of (un-physical) energy transfer 
from the DM particles to the baryon gas (Steinmetz 
& White 1997). The effect of two-body heating is ex- 
pected to decrease with the number of DM particles in 
the simulation, so the trend with resolution in all codes 
qualitatively suggests that at least part of the different 
temperature below T < 10^ K is related to this effect. 
However, the evolution of gas temperature with resolu- 
tion in ENZO runs is extremely small compared to all 
other codes. 

Interestingly, a similar trend was noticed by 
O'Shea et al. (2005), by comparing the temperature 
distributions obtained with GADGET2 and ENZO 
(both using the PPM version of the code, or its formu- 
lation with artificial viscosity, i.e. ENZO-ZEUS). The 
authors suggested that the reported trend were consis- 
tent with an increasing action of the effective viscosity 
employed in the hydro solver of the three codes, going 
from ENZO-PPM to ENZO-ZEUS to GADGET2. This 
explanation is also likely in o ur case; we will come 
to this point again in Sect. 16. 5[ in connection with the 
study of phase diagrams for the shocked cells/particles 
in the various runs. 
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Figure 9. Mass weighted profiles of gas density {left column), gas temperature {center column) and gas clumping factor {right column) for Cluster A at various 
resolutions. GADGET runs are in the upper row, TVD runs are in the middle and ENZO runs are in the bottom row. Vertical dashed lines show the minimum 
radius enclosing the minimum mass suitable for convergence studies, as introduced in Sec|4] 




Figure 10. Volume weighted profiles of gas entropy (in arbitrary code units) for Cluster A at various resolutions. The vertical dashed lines show the minimum 
radius enclosing the minimum mass suitable for convergence studies, as introduced in Sec|4] 



5.3 Properties of Galaxy Clusters 

Differently from the case of gas density and gas tem- 
perature distributions in the whole simulated volume, 
for which large statistics is available, our setup does 
not allow us to study the convergence with resolution 



of cluster statistics for a large number of objects. Given 
the minimum requirement of mass and spatial resolu- 
tions outlined in the previous Sections, we must expect 
that only a few galaxy clusters in our (100Mpc//i)^ 
box are sampled with enough particle/cells to allow 



12 EVazza, K.Dolag, D.Ryu, G.Brunetti, C.Gheller, H.Kang, C.Pfrommer 



the monitoring of thermodynamical distributions in- 
side the virial radius, being free from resolution effects, 
namely the two most massive clusters within the sam- 
ple: 

• cluster A: a system of total mass M = 1.36 • 
10^^ M^/h and i?vir = 2.32Mpc//i, in a fairly relaxed 
dynamical stage; 

• cluster B: a system of total mass M = 1.64 • 
10^^ M^/h and i?vir = 2.47Mpc//i, in an ongoing 
merger phase. 

We preliminary checked that the total masses at all 
resolutions and in all codes are in agreement within a 
^ 6 per cent level within i?vir, so that the general pa- 
rameters defining the systems are nearly identical in all 
investigated resolutions. 

However, we still have a minor source of scatter 
in the detailed comparison of data, given by the fact 
that the different codes adopt different time stepping 
criteria, and even if the cosmic time of the outputting of 
data is formally the same, tiny differences of the order 
of a few ^10 Myr can be expected in the data. This 
is expected to be problem only for the comparison of 
small scales in the cluster profiles, for which a perfect 
synchronization is impossible. 

This is issue is particularly relevant for the cluster 
merger B: at z = the exact positions of the thermal 
features linked to the merger event are spread at dif- 
ferent distances from the cluster center, as an effect of 
tiny differences in the internal timings of the codes. In 
the case of the relaxed cluster A, the spatial locations 
of sub-clumps is much more similar in all codes. 

In Fig. [9] (upper panels) we show the mass- 
weighted profiles of mass density, temperature and of 
the gas mass clumping factor, Sp =< p'^ > / < p , 
for cluster A at z = 0. 

We define here the mass density profile as 

rrii/Vsheih where m[ is the mass associated to each 
particle/cell in the simulation, and l^heii is the volume 
of each radial shell along the radius. The profile thus 
defined is independent of the differences in the proper- 
ties of clumping within each shell, and allows us to in- 
vestigate how the matter is distributed in the the differ- 
ent simulations. The computation of the clumping fac- 
tor then provides the complementary information about 
the distribution function of gas matter within each ra- 
dial shelE 

In this case, the weighting by gas mass ensures that 
the profiles are closely related to the thermal energy of 
the gas inside clusters, which in turn depends on the 



^ We notice that constructing the radial profiles at large distance from the 
center of clusters can be affected by tessellation problems in the case of 
SPH runs, if the smoothing length of the particles is large compared to the 
width of the shell used to compute the profile. The discreteness of grid cells 
(whose edges may intersect more than a single radial shell) may be regarded 
as a small source of uncertainty for the computation of the radial profiles in 
the lower resolution grid runs. Correcting for these effect is non trivial, and 
complex tessellation techniques may be adopted in order to minimize the 
above effect. We notice however than the trends reported in our work are 
generally much larger than the uncertainties associated with these issues. 



Statistics of energetic and low Mach number internal 
shocks (see also Sect.©. 

The profiles of density and temperature converge 
with resolution rather steadily, with an agreement bet- 
ter than a ^ 20 per cent between the profiles at all 
radii, when different resolution are compared. This is 
reassuring, since the combination of the above pro- 
files gives the profiles of the thermal energy distribu- 
tion within the clusters, and this is a rather well con- 
verged finding in all codes. On the other hand the pro- 
files of the gas clumping factor shows a much slower 
convergence even within each code, with sizable evo- 
lution at all radii from the cluster center. In all runs 
the clumping factor increases with radius, and reaches 
< 6p >^ 10 outside of i?vir- At the best available res- 
olutions the self-convergence between for each code is 
yet to be reached, despite the fact that the profile of gas 
matter density is much better behaved. 

In Figure [10] we report the volume weighted pro- 
files of the entropic function (S = T/p'^/^) for each 
particle/cell, for all codes and resolution. The weight- 
ing by the volume here is chosen to focus more on the 
trend of the entropy associated with the smooth, vol- 
ume filling accretion around clusters. Compared to the 
more standard "entropy" profile, based on the ratio be- 
tween temperature and p^/^ profiles, we consider the 
profile of the entropic function more useful to charac- 
terize the tiny differences of entropy which could be 
very locally associated with different dynamical accre- 
tion pattern in the different codes. 

The study of mass-weighted entropy distributions 
in a smaller re- s imu lation of this project will be dis- 
cussed in Sect. I5.4[ where we investigated the en- 
tropy generation associated with the clumps of mat- 
ter in clusters. In this case, GADGETS runs are those 
characterized by the slowest resolution compared to 
grid methods, and they also present a peak of the en- 
tropy gradient at a significant larger distance compared 
to TVD and ENZO runs. We also report the interest- 
ing trend that, compared to grid methods, the increase 
in spatial resolution causes a significant smoothing of 
the entropy jump in GADGETS clusters (see Fig. flTI) . 
At the best available resolution, the full-width-half- 
maximum of the entropy "jump" in grid methods is 
significantly smaller than in GADGETS (^ 2i?vir in 
TVD and ENZO versus - 3i?vir in GADGETS). To 
check if differences in the clumping of gas matter is re- 
sponsible for the above differences, we also computed 
the profiles for the 256 GADGETS run by consider- 
ing only the 50 per cent less dense particles (Fig. [TT1) . 
but no significant differences can be found. This dif- 
ference in GADGETS can only be partially explained 
by SPH smoothing effects, since the observed broaden- 
ing is considerably larger than the smoothing length at 
these over densities. The dynamics of shock waves on 
large scale accretion pattern around clusters are how- 
ever expected to play the major role here: we will fur- 
ther explore this issue in Sect.O 

A second interesting feature of entropy profiles 
is the hint of a flattening of the entropy profile at 
^ 0.3i?vir in clusters simulated with ENZO compared 
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to GADGETS runs. This is in line with a number of 
existing results in the literature (e.g. Frenk et al. 1999; 
Wadsley et al. 2008; Tasker et al. 2008; Mitchell et al. 
2009), even if the grid resolution here is too coarse to 
show conclusive evidence. However, tests employing 
efficient adaptive mesh refinement with ENZO have 
recently shown that that the extreme flatness of the en- 
tropy profile in these cluster runs inside 0.1i?vir is a 
very robust feature against numerical and mass/spatial 
resolution effects (Vazza 201 1). 

Based on the literature, it seems likely that the dif- 
ferences in the inner entropy profiles are produced by 
the different integrated mixing role played by artificial 
viscosity which is enhanced in grid codes compared to 
SPH (e.g. Wadsley et al. 2008; Mitchell et al. 2009). 
With our setup we tested in detail the way in which 
entropy is advected inside clusters in ENZO and GAD- 
GET with a furt her r e-simulation, discussed in the Sec- 
tion below (SecEl 

On the other hand, it is very likely that the leading 
mechanism which sets the shape of the entropy distri- 
bution beyond i?vir is the action of shock waves. Grav- 
itationally induced motions of gas matter are the lead- 
ing drivers of shock waves in these simulations, and 
therefore a detailed analysis of the distribution of gas 
matter in the outer shells of simulated clusters is help- 
ful to understand the reported differences. While we 
defer to Sect. [6l a detailed study of the morphologies 
and statistics of accretion shock around clusters, while 
here we study in detail the simple gas matter distribu- 
tion in the outer cluster regions, comparing different 
codes and resolutions. 

The panels in Fig. [121 present the mass-weighted 
and volume-weighted distribution of gas matter within 
the radial shell 1.5i?vir ^ ^ ^ 2i?vir outside of cluster 
A, for the same runs of previous Figures. 

As expected, the volume weighted distributions 
show that the grid codes are able to resolve more struc- 
tures (e.g. smooth filaments of gas) in the low den- 
sity regions; on the other hand it can be seen in the 
mass weighted distributions that GADGETS resolves 
much more collapsed objects, which are absent in the 
grid codes. This corresponds to the larger number of 
gas clumps that can be visually seen in the projected 
maps of Fig. O However, such material in grid codes 
produces also an excess of baryon gas in the range 
10~^^ gjcrv? ^ p ^ 10"^'' gjcrv? , compared to 
GADGETS. These two excesses in grid codes and in 
SPH produce signals of a similar order, which explains 
why the average clumping factors reported in Fig. [9] 
are quite similar, despite the fact that the differential 
distribution have rather different shapes. The differen- 
tial distribution of gas matter in the outer shells pro- 
vide a preliminary suggestion that the shock waves as- 
sociated with these accretions can be significantly dif- 
ferent in the two methods. Indeed, a larger contribu- 
tion from stronger shocks (driven by "smooth", rather 
than "clumpy" accretions) should be expected on av- 
erage in grid codes, at the same radius. The larger en- 
tropy jumps associated with these strong shocks around 
smooth accretions may then well explain the differ- 
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Figure 11. Cross comparison of the volume-weighted gas entropy profiles 
(in arbitrary code units) for cluster A {Jeft column) and for cluster B {right 
column). GADGETS runs at 128 are reported in dot-dash, while the 256 
runs are in solid; the long dashed lines report the profiles for GADGETS 
runs at 256, but considering only the 50 per cent less dense particles. 



ences of shape in the outer entropy profiles of cluster 
A and cluster B. For recent works employing higher 
spatial and mass resolution to characterize in detail the 
clumping and azimuthal scatter properties of gas mat- 
ter in the outer region of galaxy clusters, we address the 
reader to Roncarelli et al. (2006), Bums et al. (2010), 
Vazza et al. (2011) and Nagai et al. (2011). The issue 
of matter clumping in the outskirts of galaxy clusters 
has also recently become a topic available to X-ray ob- 
servations (e.g. Simionescu et al. 2011; Urban et al. 
201 1), and therefore the predictions of different numer- 
ical methods, even at ^ i?vir, are going to be likely 
tested with observations in the next future. 



5.4 A Test with Tracers 

In order to analyze in a more conclusive way the differ- 
ences in the entropy profiles of cluster, we performed 
a re-simulation study which followed in detail how the 
entropy of gas is build inside one massive clusters dur- 
ing its evolution. 

To this goal we simulated a smaller volume of side 
40Mpc//i, whose initial conditions were produced in 
a similar way as in Sec|3l in this case an even larger 
normalization for the matter power spectrum parameter 
were used, ag = 1.6, in order to form a M ^ IO^^Mq 
cluster inside this small volume0. 

Since the entropy profiles of grid codes were found 
to be very similar, for simplicity we tested here only 
the and ENZO run with 256^ cells (corresponding to a 
spatial resolution of 156kpc//i with a GADGETS run 
with 256^ gas particles. 

We are interested in the evolution of gas entropy 
linked to the matter accretion history of the cluster, 
and we identified all gas sub-halos in place at z = 1 
outside of the main cluster in the volume, and we fol- 
lowed their evolution in time. The location of their cen- 
ters (based on a spherical over-density halo finder) is 



^ The initial conditions for the 40Mpc//z box, at different resolutions, can 
be found at this URL: |http://canopus.cnu.ac.kr/shocks/casel/ 1 
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Figure 12. Mass -weighted (solid lines) and volume-weighted (dashed) distributions of gas density within the shell 1.5 ^ r/R^ij- ^ 2 around cluster A, for 
all simulated runs. 




Figure 13. First two columns: projected map of gas density (colors) and tracers positions for 4 time steps in ENZO 256 run, for the test simulation described 
in Sec l5.4| Last two columns: same as in the first 4 panels, but for GADGETS 256 run. The side of the images and the line of sight are comoving 12Mpc//i 
in all cases. 



in agreement in both simulations within a 200kpc//i 
accuracy. We selected all particles belonging to the 3 
sub-halos in GADGETS runs, while (mass-less) tracer 
particles were placed inside the corresponding cells 
in ENZO run. The distribution of tracers was gener- 
ated using an number density profile corresponding to 
a King profile, using a sampling of ^ 0.1 of the cell 
size. We checked that the final tracers distributions are 
statistical independent of the particular profile adopted 
for the initial generation (see also Vazza, Gheller & 
Brunetti 2010). 

The gas tracers in ENZO were then evolved by 
updating their positions according to the underlying 
Eulerian velocity field, with the same procedure of 



Vazza, Gheller & Brunetti (2010). In summary, the 
three-dimensional velocity field was interpolated at the 
location of tracers using a Cloud In Cell kernel, and 
the positions were updated every 2 time steps of the 
simulation with a first-order integration. The entropy 
assigned to the tracers at each time step corresponds 
to the entropy of the cells where each tracer sits at the 
time of observation. 

The visual inspection of projected tracers/SPH par- 
ticles positions as a function of redshifts (Fig. fT3l) 
clearly shows that the accretion of gas clumps is a dif- 
ferent process in the two runs. 

Even if the initial positions of the clumps centres 
are equal down to the cell resolution, soon after their 



A Comparison of Cosmological Codes 15 







Figure 14. Radial profiles of the physical gas entropy (in arbitrary code units) corresponding to all panels in Fig. [13] The soHd lines show the average mass 
weighted entropy profiles for the complete GADGETS/particles (left 4 panels) and ENZO/cells distributions (right 4 panels), while the overlaid colours show 
the contribute from particles/tracers initially located within the 3 selected sub-halos. 



accretion through i?vir their trajectories differ consid- 
erably: the particles from sub-halos in GADGETS soon 
mix with the main cluster atmosphere after accretion, 
and most of the particles from sub-halos end up in the 
dense and low-entropy cluster core. In ENZO the trac- 
ers mix more slowly at the beginning, and most of ac- 
creted gas component is bound to the infalling clumps 
even after the crossing i?vir- In particular, most of trac- 
ers initially located in two clumps (colored in blue and 
in gray) never penetrate inside the core of the main 
cluster, but find their selves settling at larger cluster 
radii, - 0.2 - 0.3i?vir. 

The analysis of the entropy profiles of the main 
cluster and of SPH particles/tracers is presented in 
Fig. [in and confirm the difference in the accretion his- 
tory of the two methods. In this case, since we are in- 
terested in the evolution of gas clumps, the weighting 
by gas density of the entropic function is adopted here. 

In GADGETS, only a fraction of the matter from 
clumps is shock heated to higher entropy, and the un- 
shocked low entropy material can be delivered to the 
low entropy center of the main cluster, where it remains 
until the end of the simulation. Already at z = 0.25 

5 Gyr after their accretion inside i?vir) the entropy 
of SPH particles from sub-halos is nearly identical to 
the entropy of the main cluster. On the other hand in 
ENZO run the gas from clumps is soon shock heated 
to higher entropy values (compared to particles in sub- 
halos in GADGETS), and it retains its entropy for a 
larger time, placing on average on radii external to the 
cluster core. In the ENZO run, there is still a relevant 
scatter in the entropy of tracers at z = 0. 1, compared to 
the main profile of the cluster, which is very different 
from GADGETS results. 

Our results suggest that the following different 



mechanisms are at work in the two methods: a) in SPH, 
accreted clumps soon loose their gas because of the in- 
teraction with the ICM of the main cluster, the entropy 
of their gas gets quickly in an equilibrium with the at- 
mosphere of the host cluster and many particle from 
the sub-halos can end up within the low entropy core 
of the main cluster; b) in PPM, accreted clumps are ef- 
ficiently shock heated while entering the atmosphere 
of the main cluster, they reach more slowly an equilib- 
rium with the average entropy of the main cluster at- 
mosphere and most of the accreted material sets to an 
higher adiabat in the cluster profile (compared to the 
SPH run), avoid to concentrate within the cluster core. 

In both cases, we observe that the shock heating 
and mixing motions following the matter accretions 
from small satellites (i.e. minor mergers) are not effi- 
cient processes in changing the overall shape of the en- 
tropy profile within the main cluster, which is already 
in place at z ^ 1. 

On the other hand, we can speculate that the differ- 
ent trajectories and thermodynamical evolution of the 
gas matter accreted by sub-clumps in the two methods 
highlights the sizable differences of transport phenom- 
ena in the two schemes, which are relevant to many as- 
trophysical topics in galaxy clusters (e.g. metal enrich- 
ment, cosmic ray transport, non-thermal emissions). 

Since we do not make use of adaptive mesh refine- 
ment in ENZO simulations here, the spatial resolution 
is too poor to study fluid instabilities and cluster turbu- 
lence (for studies of tracers in high resolution ENZO 
runs with adaptive mesh refinement, see Vazza, Gheller 
& Brunetti 2010 and Vazza 201 1). However we notice 
that at this point it is clear that the flatness of inner 
cluster entropy profile generally found in PPM codes 
is not a product of employing AMR itself, but it is a 
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more fundamental feature linked to shocks and mixing 
inside clusters. 

In their seminal work Mitchell et al. (2009) investi- 
gated the production of cluster entropy in a binary clus- 
ters merger with GADGET and the PPM code FLASH 
(Fryxell et al. 1998), and found that the most important 
factor which produces the differences seen in the two 
numerical methods is the early mixing of entropy dur- 
ing the collision of cluster cores, driven by fluid insta- 
bilities, which is much more pronounced in PPM than 
in SPH. Our test here shows that the way in which fluid 
instabilities and shocks follow the accretion of smaller 
subunits of cluster also differ in the two approaches, 
and lead to dissimilar entropy tracks for the accreted 
gas. 



6 SHOCK WAVES IN COSMOLOGICAL SIMULATIONS 

Many of the differences previously found between the 
codes, such as the temperature structures in low den- 
sity environments and entropy distributions in the in- 
nermost and in the outer regions of clusters, are likely 
connected to the dynamics of matter accretion pro- 
cesses in the accretion regions of large scale structures. 
In these regions, the activity of strong shock waves is 
the leading driver of thermalization, entropy genera- 
tion and possibly of cosmic ray acceleration in large 
scale structures, (e.g. Ryu et al. 2003), via the diffu- 
sive shock acceleration mechanism (e.g. Blandford & 
Ostriker 1978). 

The numerical modeling of shock waves is among 
the most important tasks that cosmological codes must 
correctly perform in run-time; several different numer- 
ical techniques, involving the use of ad-hoc numerical 
viscosity (as in SPH) or the solution of the Riemann 
problem through explicit methods (as in PPM or TVD), 
have been adopted for this task (e.g. Dolag et al. 2008 
for a review). 

All these methods generally perform well in the 
case of rather simple shock problem (e.g. Tasker et al. 
2008), while their performances in the very complex 
environment of large scale structure simulations are 
more uncertain. To date, no detailed comparison of the 
statistics of shocks developed in the various numerical 
method have ever been published; our sample of runs 
thus offers the optimal framework to test the outcomes 
of the different methods in the cosmic volume. 

In order to readily compare the statistics of shocks 
in each simulations, a shock finding method is needed 
to detect and measure the strength of shocks in the sim- 
ulations. To this end we start by presenting the shock 
detecting method explicitly developed to work on each 
specific code in our project. 



6.1 Shocks Capturing Algorithm 

The Rankine-Hugoniot jump conditions allows one to 
evaluate the shock Mach number, M, from the ther- 
modynamical state of the pre-shock and post-shock re- 
gions (under the assumption of a pre-shock medium at 
rest and in thermal and pressure equilibrium). If the 



adiabatic index is set to 7 = 5/3 one has the well 
known relations (e.g. Landau & Lifshitz 1966): 
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with indices 1, 2 referring to pre and post-shock quan- 
tities, respectively, and where the entropy S is S = 

In practice measuring M of shocks in cos- 
mological simulations is more problematic than in 
this ideal case: matter falling in the potential wells 
drives chaotic motions and the temperature distribution 
around shocks is usually patchy due to the continuous 
accretion of cold clumps and filaments into hot halos. 
These complex behaviors establish complex pattern of 
pre-shocks velocity, temperature and density fluctua- 
tions which makes problematic to measure Ranking- 
Hugoniot jumps in a clean way. To overcome this prob- 
lem, detailed analysis strategies have been conceived 
over the last years, with the goal of recovering the 
measure of M in fully cosmological simulations in the 
most accurate way. 



6.1.1 The Temperature Jumps Method - TJ 

The analysis of jumps in temperature is a powerful way 
of measuring the strength of shocks in Eulerian cosmo- 
logical simulations, and its application was first dis- 
cussed in Miniati et al. (2001), with a more sophisti- 
cated formulation in Ryu et al. (2003). The cells host- 
ing a possible shock pattern are preliminarily tagged 
by two conditions: 

• VT • V5 > 0; 

• V • V < 0. 

The additional condition on the strength of the 
temperature gradient across cells is also customary re- 
quested: 

• I MogT |> 0.11; 

(specifically | AlogT |> 0.11 filters out shocks with a 
Mach number M < 1.3, Ryu et al. 2003). 

It is customary to simplify the process of identi- 
fication of shocked cells by using a one-dimensional 
procedure applied successively in three orthogonal di- 
rections. In the case of multiple shocked cells in close 
contact, the center of shocks, which can be spread 
across 2-3 zones, is placed where V • v is minimum. 
Then the Mach number is calculated based on Eq 121 
where T2 and Ti are the post and pre-shock temper- 
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Figure 15. Projected maps of shocks (in logM) for a slice of 75Mpc//i in the simulated volume for the most resolved runs of the sample (left:ENZO at 512^ ; 
center: GADGET at 256^; right: TVD at 512^). We adopt a weighting by volume for each particle/cells, and a fixed width of ~ 550kpc along the line of 
sight in all maps. 
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Figure 16. Map of shocks (in logM) for a slice with the side of 25 Mpc/h through the center of cluster B. The first row reports the results of the TVD runs 
and the Temperature Jump shock finder as a function of resolution, the second row reports the results for the PPM runs and the Velocity Jumps shock finder, 
the third row reports the results fro the SPH runs and the Entropy Jump shock finder. From the left to right column, the width along the line of sight is 2200, 
1100, 550 and 275kpc respectively. 
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Figure 17. Volume weighted number distributions of shocks at all resolutions and for all simulations. In the right panel also the results from the grid codes (at 
the maximum available resolution) are shown for comparison 



ature across the shock regiorQ. In the following Sec- 
tions, we will refer to this method as to the TJ method. 

In this work, we applied the TJ method following 
the original formulation of Ryu et al. (2003), with the 
exception that we do not employ the temperature floor 
of To = 10^ customarily used to mimic the effect of 
re-ionization, in order to readily compare with the out- 
comes of the other simulations of the project. 



6.1.2 The Velocity Jumps Method - VJ 

A similar approach, based on the post-processing anal- 
ysis of velocity jumps across cells in grid simulations 
was proposed in Vazza, Brunetti & Gheller (2009) for 
the analysis of ENZO simulations. Conservation of 
momentum in the reference frame of the shock yields: 

p^vx = 92^2, (4) 

with the same notation used in EqsiT]-[3l In the ideal 
case in which the pre-shocked medium is at rest and 
in thermal and pressure equilibrium, the passage of a 
shock with velocity Vs leaves a /\v in-print as a ve- 
locity difference between the shocked and pre-shocked 
cells. In the lab frame a relation holds between l\v and 
M, which can be obtained by combining Eqn. |4] with 
Eqn.n 



3 1 -M^ 



(5) 



where Vs = Mcs and is the sound velocity com- 
puted in the pre-shocked cell. 

The procedure to identify shocks in 3-D with the 
VJ method follows these steps: 

• candidate shocked cells are selected as V • v < 
(calculated as 3-dimensional velocity divergence); 



^ We note that Skillman et al. (2008) pointed out that the application of a 
split coordinate approach to the TJ method may lead to an overestimate in 
the number of shocks, compared to an unsplit TJ method, in ENZO AMR 
simulations. The bulk of the thermalized energy at shocks, however, is only 
marginally affected by the above differences 



• if more candidate shocked cells are found to- 
gether, the one with the minimum V • v is considered 
as the shock center; 

• the three Cartesian axes are scanned with 1-D 
sweeps and Avx,y,z jumps along the axis of scan are 
measured, between cells located at a A/ distance on 
opposite side of the shock center. In ENZO PPM we 
can safely use A/ = 1, therefore M is measured across 
3 cells (e.g. Vazza, Brunetti & Gheller 2009 for a de- 
tailed discussion). 

• the sound speed is taken from the cell in the 
tagged patch which shows the lower temperature, and 
based on this the Mach number along each direction is 
computed from Eqn. [51 

• we finally reconstruct the 3-D Mach number in the 
shocked cell with M = {Ml ^ Ml ^ Mlfl^. 

In the following we refer to this procedure as the 
velocity jump (VJ) method. 

Vazza, Brunetti & Gheller (2009) reported overall 
consistency between VJ and TJ method in ENZO sim- 
ulations with fixed grid resolution, with minor differ- 
ences in the most rarefied environments. In Vazza et al. 
(2009) and Vazza et al. (2010) the application of the VJ 
method is extended to ENZO runs with Adaptive Mesh 
Refinement. 

The application of a qualitatively similar method, 
working on the velocity field of SPH particles in GAD- 
GET3 simulations, has also been presented by Hoeft et 
al. (2008) 



6.1.3 The Entropy Jumps Method - EJ 

A method to measure the Mach number of gas flows 
in GADGET runs was presented in Pfrommer et al. 
(2006). In this method, a run-time algorithm monitors 
in run-time the evolution of entropy for each particles, 
and from the entropy jump (in time) the Mach number 
of the shock can be inferred. 

The instantaneous injection rate of the entropic 
function due to shocks for each SPH particle is 
dA(S)/dt, where A is the entropic function, defined 
by P = A{S)p^ (where P is the gas pressure). If 
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the shock is broadened over a scale of order the SPH 
smoothing length f^h (fh ^ 2 is a factor which has to 
be calibrated with shock-tube tests), one can roughly 
estimate the time it takes the particle to pass through 
the broadened shock front as At = //^/i/ v, where v can 
be approximated with the pre-shock velocity vi. As- 
suming that the present particle temperature is a good 
approximation for the pre-shock temperature, it is pos- 
sible to replace vi with Mici. 

Based on these assumptions and using AAi 2^ 
AtdAi I dt, the jump of the entropic function of the par- 
ticle crossing a shock will be: 
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where, using Equation [T] and |2] one has: 
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that combined with Equations [6] and |71 
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The right-hand side of Eqn. [9]can be estimated in- 
dividually for each particle, and Eqn. [9] allows to esti- 
mate their Mach number (see Pfrommer et al. 2006 for 
details). In the following we will refer to this method 
as EJ method. 

The EJ method has been applied in a series of pa- 
pers to characterize shocks on the fly, inject CRs with 
a Mach number-dependent acceleration efficiency, ac- 
count for the non-linear back reaction of the CR pres- 
sure on the hydrodynamics and following the trans- 
port of CRs during GADGETS simulations of cosmo- 
logical structure formation, galaxy and galaxy cluster 
formation (Pfrommer et al. 2006, 2007, 2008; Pfrom- 
mer (2008); Jubelgas et al. (2008); Pinzke & Pfrommer 
2010; Pinzke et al. 2011). In our work here, the origi- 
nal EJ scheme has been applied in run-time to GAD- 
GETS runs, and the measured distributions of Mach 
numbers for the gas particles have been analyzed in 
post-processing. 



6.2 Shocks Maps and Morphologies 

We measured the strength of shocks in our simulations, 
by applying the TJ method in post-processing to TVD 
runs, the VJ method in post-processing to ENZO run 
and the EJ method in run time for GADGETS runs. 

The panels in Figure [15] show the large scale pat- 
tern of shock waves for a thin slice (of 550 kpc) in the 
simulated box at z=0, for the best available resolutions 
in all codes. Only for display purposes, the Mach num- 
bers measured in GADGETS have been interpolated 
onto regular grids with resolution corresponding to a 
256^ mesh. 



Even at the best available resolution, the mor- 
phological distributions of shocks in the various run 
looks less similar than what is generally found for the 
density-weighted maps of temperature (Fig.|4ll5l). In all 
runs innermost region of clusters and filaments hosts 
only weak shocks, M ^ 2 — 5, while the strongest 
shocks are located outside cosmic structures. However, 
the strong external shocks are very sharp and regular in 
grid codes, while they seem to be grouped in clumps 
in GADGETS. While in GADGETS runs the shocked 
structures are rather volume filling (due to the smooth- 
ing kernel in less dense regions), in both grid methods 
the shocks outside clusters are regular surfaces with ra- 
dius of curvature ^ 3 — lOMpc, with a very small vol- 
ume filling factor. 

We notice that this difference between SPH and 
grid methods depends on the different resolutions out- 
side i?vir, however the general trend is that when the 
spatial and mass resolution of DM particles is in- 
creased, the differences between grid codes and SPH 
are even more sizable. 

This is shown in the panels of Fig. [161 which zoom 
into the cluster region at the center of the cosmological 
box. Looking at the strong external shocks in the upper 
left sector of the cluster, one can see that these features 
become increasingly sharper and more regular in grid 
methods, while they become stronger and more clumpy 
in GADGETS runs. On the other hand, the trend with 
resolution inside of the cluster is quite similar in all 
codes, with increasingly thinner and weaker shocks as 
the resolution is increased. 



6.3 Mach Number Distributions 

The volume distribution of Mach numbers in the cos- 
mological volume is a simple statistical proxy that al- 
lows us to readily compare the different shock finder 
and underlying simulations. However, they cannot be 
directly translated into observational quantities, and 
therefore their study is just intended to be a useful to 
cross-check of numerical implementations, rather than 
a physical test. 

Figure [TT] shows the volume-weighted distribu- 
tion of shocks Mach number from all runs using our 
projects. 

At the best available resolution, the distributions 
from the different methods are quite similar, showing 
a peak of shocks at M ^ 1.5 and a steep decrease at 
stronger shocks. Compared to the peak, the average fre- 
quency of M > 1000 shocks is - 10"^ in GADGET, 
and - 10-^ in ENZO and TVD. 

GADGET runs present the best degree of self- 
convergence, with very little evolution between runs 
64^ and 256^. The VJ methods applied to ENZO runs 
on the other hand shows the slowest degree of evolu- 
tion, with a particularly poor performance at the 64^ 
run; this is due to the difficulty of removing baryon 
bulk flows from velocity jumps associate with shocks 
at very low grid resolutions. The TJ method present a 
noticeable self-convergence at all resolutions, although 
a the 64^-128^ run present a different convexity in the 




range 10 ^ M ^ 100 (where the contribution from 
internal and external shocks takes place), similar to the 
converged findings of the EJ method applied to GAD- 
GETS. 

In both grid codes, the increase of resolution al- 
ways cause a progressive weakening of the strongest 
shocks in the most rarefied environments; also the 
bump of external shocks is progressively shifted to- 
wards lower M. 

We notice that at the best available resolution here, 
the convergence all simulations (and most significantly 
in grid codes) is not yet reached, even if it looks ap- 
proaching; the same is true also for the distribution of 
thermal energy flux across shocks (Sect l6.4l) . Based on 
the tests in the literature, run with these same codes 
(e.g. Ryu et al. 2003; Skillman et al. 2008; Vazza, 
Brunetti & Gheller 2009; Vazza et al. 2009) one can see 
that a very good convergence (i.e. better than a ^ 10 
per cent level) in the most important shock statistics is 
expected to for a spatial resolutions of ^ 50 — lOOkpc, 
which are below our best resolution here. However the 
trend with resolution is usually very regular, and the 
differences reported here are significant, despite the 
fact that a small evolution with resolution may still be 
present. We also remark that an additional and unavoid- 
able source of difference with resolution is due to the 
ways in which the shock-finder methods work, because 
that the dependence on resolution of the different ther- 
modynamical jumps used for the computation can be 
different, especially for very coarse resolution. 

We also notice here that the modeling of a re- 
heating UV radiation from from massive stars and 
AGNs is crucial for a realistic estimate of the baryon 
gas temperature outside of cosmic structures (e.g. 
Haardt & Madau 1999). In order to measure realistic 
Mach number in the rarefied universe outside clusters, 
groups and filaments a re-ionization temperature back- 
ground is usually accounted in simulations, either in 
post-processing (Ryu et al. 2003; Skillman et al. 2008; 
Vazza, Brunetti & Gheller 2009) or in run-time (Pfrom- 
mer et al. 2006; Vazza et al. 2010). In this case the 
minimum temperature in all si mula tions is set by the 
low temperature floor (see Sect. 15.21) : however the dif- 
ferences in the values adopted in the different codes 



(from 1 K in the case of ENZO, to 24 K in the case of 
GADGET3) cannot account for the sizable differences 
in the distribution of Mach numbers. 

The differences between the methods are high- 
lighted when we plot the volume weighted average 

Mach number of shocks, M, as a function of gas den- 
sity (Fig. [18]). The results of the different codes are 
consistent only for p/pcr ^ 10 regions (typical of 
the outskirts of galaxy clusters and filaments), with 

M ^ 2. At lower densities we report the following 

trend: in SPH M is smoothly increasing moving to- 
wards lower density regions, while in grid codes the 

transition of M moving to lower densities is very sharp, 

and causes a net increase of M by 2 orders of mag- 
nitude in both grid methods. These large differences 
in the range p/pcr < 10 mirror the different thermal 
structures of baryons in the outermost regions of LSS 
in grid codes and in SPH (SecH]). In these environ- 
ments, the self convergence in grid codes is not yet 
reached even at the best available resolution 195kpc//i. 



6.4 Energy Distributions. 

The thermal energy flux across each shock in the sim- 
ulations is measured as: 

/th = S{M) . (10) 

where ppre is the pre-shock gas density and S{M) 
is a monotonically increasing function of M which fol- 
lows from Rankine-Hugoniot jump conditions, whose 
formula can be found for instance in Kang et al. (2007). 

In the TJ and in the VJ methods this quantity is 
computed in post-processing based on the shock direc- 
tion, while in the EJ method /th it is measured in run- 
time. 

We remark that in all 3 methods, the numerical 
recipes to compute the effective thermalization at the 
post-shock are tuned to remove the effect of adiabatic 
compression of the gas in the post-shock region, which 
can provide sizable additional thermalization in the 
regime of weak shocks (see Ryu et al. 2003; Pfrommer 
et al. 2006). 
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Figure 19. Maps of projected of the thermalized energy flux at shock waves (in units of log[erg/s]) across the whole simulated volume (top panels) and for a 
sub-region of side 25 Mpc/h centered on the most massive galaxy cluster of the sample (bottom panels). 



We report in Fig. \T9\ the projected map of /th 
across the simulated volume in the best resolved runs 
(top panels), and for a zoomed region of 25Mpc//i 
(bottom panels). We also report for comparison with 
the SPH run the corresponding ENZO 256^ run. The 
flux coming from the innermost cluster region looks 
morphologically similar in all cases, with a compact 
and spherical "envelope" of energetic shocks concen- 
trated inside the virial volume of halos. The differences 
are more sizable at the scale of filaments and in the 
outer region of clusters, where we notice very sharp 
shock surfaces even in projection in grid methods, 
while much smoother pattern are found in GADGETS, 
with external accretion shock extending at larger dis- 
tances from the center of clusters. This effect mirrors 
the corresponding distribution of gas entropy at larger 
scales, which we reported in the analysis of the radial 
profile of the entropic function, in Sec {5.3[ The zoomed 
images of Fig. [19] additionally shows that complex in- 
tersections of merger shocks are modeled inside the 
over-dense regions in grid codes, while very smooth 
distribution appears in the projected GADGET maps. 
Taking as a reference the ENZO run with 256^, we see 
that the above differences are not trivially due to reso- 
lution effects, since the large scale shock patterns in the 
grid code do not significantly get smoother or shift in 
position even if the resolution of the simulation is made 
coarser. The differential distributions of /th for all runs 
is reported in Fig.[20l In this case the contribution com- 
ing from the low density regions is fairly negligible and 
results are found to be in an overall good agreement. 
As in the case of number distributions, the EJ method 



presents the largest degree of self-convergence, and the 
VJ presents the slowest degree of self-convergence. 

The grid codes present the clear trend of process- 
ing less thermal flux at M >> 10 shocks when reso- 
lution is increased, while in SPH slightly more energy 
flux is processed at strong shocks when resolution is in- 
creased (although this amount is negligible compared 
to the peak of thermalization in the box). In the bottom 
panels of the same Figure, we also show the cumulative 
distribution for the same run, normalized to the total 
flux inside the cosmic volume for each run. 

At their best resolution, all codes agree in sev- 
eral important findings: a) the peak of thermalization 
is found at M ^ 2, consistent with most of previous 
works in the literature (e.g. Ryu et al. 2003; Pfrom- 
mer et al. 2006,2007; Vazza et al. 2009,2010; Skill- 
man et al. 2008); b) the general shape of the distribu- 
tions is quite similar, with a steep power-law behaviour, 
dlogfth/dlogM ^ . The slope is a ^ 3 in grid 

codes and a ^ 2.5 in GADGETS runs; this is steeper 
compared to the findings in the literature, because we 
are not modelling here the re-ionization background, 
c) the cumulative distributions for M < 10 shocks are 
very similar in all codes, and only ^ 1 per cent of the 
total thermal flux inside the cosmic volume belongs to 
shocks with M > 10. These findings suggest that, de- 
spite sizable differences in the shapes and statistics of 
strong external shocks in the accretion regions of large 
scale structures, the bulk of the energetic properties of 
shocks within the cosmic volume is a rather well con- 
verged answer from cosmological simulations. 
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Figure 20. Top panels: differential distributions of the thermalized energy flux through shocks at all resolution and for all codes. Bottom panels: cumulative 
distributions for the same runs. 



6.5 Phase Diagrams for Shocked Regions. 

To pinpoint the differences between the codes, we find 
it useful to extract the phase diagram of shocked cells 
for the various runs within the total cosmic volume. 
Panels in Fig. [2T] and [22l show the flux-weighted mean 

Mach number, M, and thermal flux (normalized to 
the total thermal flux in the cosmic volume) for the 
shocked cefls of runs 64^, 128^ and 256^. 

In grid codes, as soon as the spatial resolution is 
large enough to model the innermost region of col- 
lapsed halos, a compact "group" of cells at M ^ 10 
is formed in the upper right comer of the phase dia- 
gram, while a much broader region of strong shocks is 
found at lower densities at across a wide range of tem- 
peratures. In GADGETS, a similar "group" of points 
corresponding to halos is formed, but it has less sharp 
contours and it smoothly extends to lower densities, 
where strong outer shocks from a concentration which 
is much narrower compared to grid codes. 

If the dissipated energy flux is concerned (Fig. [22b. 
again less disagreement is found among codes. At all 
resolutions, about the ^ 90 per cent of total dissipated 



energy in the box is found at cells with p/ per ^ 10^ 
and T > lO^K. 



One should expect a high degree of convergence 
in the statistics and morphologies of energy dissipating 
structures in the three codes: indeed the main sources 
of heating in these adiabatic runs are shocks, and the 
cro ss co mparisons in the previous Sections (Sec i5.ll - 
Sec i5.3l) have shown that most of the thermal properties 
of halos are in good agreement. 

On the other hand, shocks are also the main source 
of entropy generation in these simulations, and we 
showed that the halos in the different codes present 
noticeable differences b oth in th e inner and outer en- 
tropy distributions (Sec i5.3ti54l) . are likely related to 
details of shocks dynamics away from the most dissi- 
pative structures in simulations. 

Fig.[23]shows the illustrative case of the scatter plot 

for the post shock entropy versus M diagram. We re- 
strict to T > 100 K regions in order to avoid any arti- 
facts due to different low te mpe rature floors adopted in 
the various codes (see Sect i5.2l 

A concentration of high entropy and weak shocks 
(in red color, in the Figure) is common to all simulated 
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Figure 21. Phase diagrams for shocked cells in the simulations, color coding shows the flux- weighted average Mach number. Additional isocontours with a 
coarse binning in M space are shown for clarity. 



data, and marks the shock energy dissipation in inner- 
most region of galaxy clusters. 

However, in grid codes a concentration of points 
is also present for M > 10^, as diagonal stripe in the 

plane (M,S). The points in this region (in blue color) 
trace external shocks, for which the post shock entropy 

is tightly correlated with M (Eqn. [3]) for strong M > 

10 shocks, leading to a 5^2 oc M^. 

This "phase" of shocked gas is almost completely 
missing in SPH runs. 

We verified that in the grid codes, the strong shocks 

following the 5^2 oc correlation are systematically 
located at the outskirts of galaxy clusters and filaments, 
while the concentration at M < 10 shocks comes from 



cells within collapsed halos. In this second case, ener- 
getic and weak shocks are unable to change the post- 
shock entropy in a relevant way, and no strong relation 

is found between S and M. Therefore, in GADGET 
entropy is released in the simulation at the same loca- 
tion of the most dissipative structures in the universe, 
whereas in both grid codes a sizable amount of entropy 
is also released at outer accretion shocks, which are not 
responsible for sizable energy dissipation. 

This suggests the important point that, although 
thermalized energy is processed in the various codes in 
a rather consistent way, the gas entropy in grid codes 
and in SPH is increased in shock structures with rather 
different morphologies and thermodynamical proper- 
ties. Considering that the production of entropy at outer 
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shocks is also respons ible for the innermost entropy 
profile in clusters (Sec i5.4L we suggest that this find- 
ings is also relevant to understand the detailed prop- 
erties of advection of matter (and possibly CR) inside 
galaxy clusters, over cosmic time. 

One possibility is that the absence of strong en- 
tropy generation at outer shocks in GADGETS is due 
to pre-shock entropy generation by to artificial viscos- 
ity (e.g. O'Shea et al. 2005), which would also be con- 
sistent with the trend reported in the temperature dis- 
tributions of Sec i5.2[ An additional effect here is likely 
the smearing of shocks at low densities in SPH, which 
makes difficult to the shock solver in GADGETS to up- 
date the particles entropy in a fully consistent way, if 
several smeared shocks merge together in the accretion 
regions. 



6.6 Shocks in Clusters and Cosmic Rays Acceleration 

Galaxy clusters are expected to be the most important 
accelerators of CR in the universe (e.g. Miniati et al. 
2001; Ryu et al. 2003; Pfrommer et al. 2007); it is 
therefore important to analyse in detail also the esti- 
mated properties of CR acceleration at shocks, in the 
most massive galaxy clusters of our simulated volume. 
In Fig. [24] we report the average radial profile of 

mean Mach number, M, for clusters A and B, show- 
ing both the results of the weighting by gas density 
(dashed lines), and by the dissipated energy flux (solid 
lines). Despite the different dynamical state of the two 

systems, we measure M ^ 2 for r < 0.5i?vir in all 
runs. Approaching the cluster virial radius, the grid 
codes show a sharp increase in the mean Mach num- 
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Figure 23. Mach versus Entropy diagrams for shocked regions of the 256^ runs. 




Figure 24. Top panels: profiles of density- weighted and energy flux- weighted average Mach number for cluster A (left panel) and cluster B {right panel). 
Bottom panels: profiles of the CR acceleration efficiency, /cR/Zth^ for cluster A {left) and cluster B {right) at the best available resolutions in all codes. 
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ber (weighted by dissipated flux), which reaches strong 

shocks, M - 10 at - 2i?vir. In GADGETS the in- 
crease in the mean shocks strength is smooth, and 

M < 3 is always found inside i?vir- The above trends 
are similar but less evident, if the weighting by gas den- 
sity is adopted. These two tre nds mirr or the trends in 
the outer entropy profiles (Sec l5.3l - [54l) . and can be ex- 
plained by noting that the medium is more clumpy in 
GADGETS runs, and that the shocks are always thin- 
ner and stronger at this location in grid codes, marking 
a very sharp the transition between large scale struc- 
tures and the rarefied Universe. 

In order to explore the possible effect played by the 
above differences in the global efficiency of clusters to 
produce the CR energy flux at shocks, we applied to 
all simulations a recipe to estimate the CR accelera- 
tion efficiency at shocks, with a standard application of 
the Diffusive Shock Acceleration theory (e.g. Kang & 
Jones 2002). According to this model, the CR acceler- 
ation at each shocks is parametrized as a function of 
the Mach number: 

fcR^v{M)-ppreM^vl/2; (11) 

where 77 (M) is a monotonically increasing func- 
tion of M, whose numerical approximation can be 
found instance in Kang et al. (2007). This prescription 
for the acceleration of CR particles is quite idealized, 
and that more recent work by the same authors take 
also into account Alfven waves drift and dissipation at 
the shock precursor (Kang et al. 2007), causing a lower 
acceleration efficiency for shocks with M < 10. Also, 
this recipe neglects the role of the re-acceleration of 
pre-existing CR, which can as well affect in a signifi- 
cant way the efficiency of acceleration at weak shocks 
(e.g. Kang & Ryu 2010). 

The bottom panels in Fig. [24l shows the radial 
profiles for the mean acceleration efficiency at shocks 
/cR//th. for cluster A and cluster B at the best avail- 
able resolutions in all codes. In the relaxed cluster A, 
the agreement is reasonably good and all codes show 
a minimum efficiency /cR//th ^ 0.1, at the cluster 
core, with a similar increasing profile up to a maximum 
of /cR//th ^ 0.7 at i?vir- Outsidc of this radius, the 
trends of grid codes and SPH largely diverge as in all 
cases reported before, and the acceleration efficiency 
in GADGETS run decreases. 

The comparison of the results for cluster B suf- 
fers of the timing issue reported in Sec l5.3[ which are 
further amplified by the non linearity of Eqn. [TT] This 
produces a large scatter from code to code in /cR//th 
inside the cluster, but approaches the same values and 
trend of cluster A for > i?vir- 

We stress that the reported differences for cluster 
B are representative of the level of intrinsic scatter that 
simulations with different numerical codes are subject 
to, which in turn adds a level of unavoidable uncer- 
tainty when estimates of CR injections from clusters 
are estimated using too small number of objects. 

It is worth stressing that the above estimate of 



CR acceleration efficiency are already at the edge, if 
not outside, of the range of permissible energy ratio 
between CR and thermal gas from gamma rays (e.g. 
Reimer 2004; Pfrommer & EnBlin 2004; Aharonian et 
al. 2008; Ackermann et al. 2010; Aleksic et al. 2010; 
Donnert et al. 2010; Pinzke et al. 201 1), radio (Brunetti 
et al. 2007; Brunetti et al. 2008) and X-ray/optical ob- 
servations (Churazov et al. 2008). Given the fairly sim- 
ple setups of the simulations considered in this work 
(e.g. no radiative processes, no re-ionization, idealized 
recipe for CR acceleration at shocks, no self-consistent 
CR feedback, coarse spatial and mass resolutions, no 
magnetic fields), this is not surprising and it suggests 
that a completely self-consistent treatment of CR, in 
presence of other important non-thermal component 
(such like magnetic fields) is needed to model obser- 
vations. 

On the other hand, these findings may also im- 
ply that the numerical implementation of the complex 
non-linear physics of non-thermal phenomena in large 
scale structures can be subject to additional uncertain- 
ties, because the basic thermo-dynamical evolution of 
accreted cosmic baryons in large scale structures is not 
yet unambiguously constrained even by rather simple 
cosmological simulations. 



7 DISCUSSION 

In this work we presented the results of a numeri- 
cal study which compares cosmological simulations 
at various resolutions, obtained with with GADGETS 
(Springel 2005), ENZO (Norman et al. 2007) and TVD 
(Ryuetal. 1993). 

The chosen simulation setup is very simple (only 
gravity forces and non-radiative hydrodynamics are 
modeled) and it is therefore particularly suitable to 
study the convergence among widely used, comple- 
mentary numerical approaches. This kind of compar- 
ison my also be helpful to explore the reasons for dif- 
ferences in the thermal and non thermal properties of 
galaxy clusters runs. 

We have analyzed in detail the properties of the 
DM distribution, thermal gas matter distribution, shock 
waves and CR acceleration efficiencies within the sim- 
ulated volume in all codes, and we highlighted all most 
convergent and least convergent findings of all codes, 
as a function of the numerical resolution and of cosmic 
environment. 



7.1 Summary of dark matter and thermal gas properties 

An overall satisfactory agreement between the 3 codes 
is found for runs with DM mass resolution better than 
^dm < 4 • IO^^Mq/Zi, in line with previous compar- 
ison works. In particular, we report a good cross con- 
vergence of the following measures: 

• the mass distribution function and baryon fraction 
for halos in the simulations are found in agreement 
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within a ^ 5 — 10 per cent, across a range of masses. 
The rate of convergence with resolution in grid codes 
is much slower than in GADGETS. These results are 
in line with the works by O'Shea et al. (2005) and by 
Heitmann et al. (2008); 

• the profiles of DM matter of halos are well con- 
verged in all codes, for all the virial volume except 
for the scales close to the gravitational softening of all 
codes, consistently with the literature (e.g. Frenk et al. 
1999); 

• the gas density distribution are agreement within 
10 — 20 per cent, for densities in the range 1 ^ p/pcr ^ 
100. High density peaks are found to be located at 
equal positions within the spatial resolution of the sim- 
ulations; 

• the gas temperature distributions are in agreement 
with a 5 — 10 per cent accuracy only for T > lO^i^ 
regions, which correspond to the typical virial temper- 
ature of the smallest halos produced in the simulations, 
in agreement with the findings reported by Kang et al. 
(1994) and O'Shea et al. (2005); 

• the gas temperature and the gas density profiles 
of the most massive clusters in the sample are simi- 
lar within a 10 — 20 cent accuracy, consistently with 
Frenk et al. (1999). Time integration of a chaotic sys- 
tem results in slightly different spatial realizations of 
substructure, in particular during mergers. This intro- 
duces an episodic source of additional uncertainty. 

On the other hand, noticeable differences are found 
in the following measures: 

• the gas density and gas temperature distributions 
for p/pcr < 1 and for T < lO^K regions are in dis- 
agreement up to 2 — 3 orders of magnitude among sim- 
ulations, even at the best available resolutions in the 
project; 

• the entropy profiles for clusters simulated with 
grid codes show a sharp peak located at ~ 2 — 3i?vir, 
while the profiles in GADGETS present a similar 
shape, but spread across a sizable lager volume; 

• the inner entropy profile of clusters simulated with 
ENZO is flat inside ~ O.li^vir, while it is steep in 
GADGETS (consistently with early results from Frenk 
et al. 1999 and more recent ones by Mitchell et al. 
2009); 

• the gas clumping within the most massive ha- 
los, and expecially in the outermost cluster regions, is 
rather different if grid codes and GADGETS are com- 
pared; 

• the time evolution of the accretion of matter 
clumps is also found to be radically different when 
ENZO and GADGETS are compared: in grid codes 
their initial entropy is substantially increased by shock 
heating, while in SPH shock heating mechanisms are 
more gentle. The accreted material is distributed at 
larger cluster radii in ENZO than in GADGETS. 



7.2 Summary of shocks properties 

Shocks were identified in all runs according to the 
shock detecting schemes specifically conceived for 
each simulation: and Entropy based method for GAD- 
GETS (Pfrommer et al. 2006), a temperature based 
method for TVD (Ryu et al. 200S) and a velocity based 
method for ENZO (Vazza, Brunetti & Gheller 2009). 

The most interesting convergent findings are: 

• the peak of thermal flux in the universe is at M ~ 
2, and originates in shocks internal to clusters; 

• the volume distribution and thermal energy flux 
distribution are very steep, and are dominated by strong 
M 100 — 1000 shocks in the external regions of large 
scale structures; 

• ^ 99 percent of the total thermal energy flux in 
the universe is processed by shocks with M < 10; 

• inside the virial radius of the most massive clus- 
ters, the density weighted profile of shocks are very 

flat, with M - 2; 

• the estimated acceleration efficiency of CR (as- 
suming Kang & Jones 2002) is small in the innermost 
cluster region, /cR//th ^ 0.1, and increases towards 
the virial radius, with /cr/ fth ^ 0.8 (however, the ab- 
solute numbers are likely to change as this recipe does 
not account for Alfven wave drift and dissipation at the 
shock precursor). 

On the other hand, the findings where we do not 
find agreement at the investigated resolutions are: 

• shocks in grid codes are morphologically simi- 
lar at all resolutions, while shocks in GADGETS show 
substantial difference at external shocks; 

• the volume-weighted mean Mach number for 
p/ Per < 10 presents different trends in each code; 

• in the vast majority of the simulated volume (out- 
side halos), shocks in grid codes show rather different 
properties in the phase diagrams (p versus T and S ver- 
sus M) compared to shocks in GADGETS. In particu- 
lar, strong accretion shocks in grid codes are associ- 
ated with large entropy jumps, while accretion shocks 
in GADGETS are not characterized by large values of 
entropy; 

• in massive clusters, grid codes produce a sharp in- 
crease of the shock strength outside i?vir, while a con- 
tinuous transition to weaker shocks is is found in GAD- 
GETS runs; 

• the CR injection efficiencies outside the virial ra- 
dius show different radial trends when grid runs and 
GADGETS runs are compared. 

7.3 Conclusions 

Overall, when cosmological numerical simulations 
with GADGETS, ENZO and cosmological TVD are 
compared within similar range of DM mass resolution, 
we report agreement better than ^ 10 per cent level 
in many statistics concerning hot, over-dense regions 
of the universe (i.e. halos, filaments). This is reassur- 
ing and it is in line with a number of previous works 
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dealing with similar topics (Kang et al. 1994; Frenk et 
al. 1999; Heitmann et al. 2008). The statistical distri- 
butions of halos masses, halos baryon fraction, density 
distributions, thermal profiles and internal shocks are 
characterized by a high rate of convergence with res- 
olution in GADGET. In the most over-dense regions, 
ENZO and TVD converge at a rather small rate, but 
produce very similar estimates at the end for most of 
investigated cases, despite the radically different hy- 
dro method they use to solve baryon gas dynamics. 
The application of Adaptive Mesh Refinement tech- 
niques is expected to further reduce the discrepancy 
between grid methods and SPH, at least in some cases 
(e.g. O'Shea et al. 2005; Tasker et al. 2008; Robertson 
et al. 2010). In the case of lower density regions (i.e. 
outer accretion regions of clusters, voids) the temper- 
ature distributions, entropy distributions, shock mor- 
phologies and Mach number distributions converge to 
rather different estimates when SPH and grid codes are 
compared. The role played by the effective viscosity 
and diffusivity of each method away from shocks may 
be partially responsible for the above differences. 

One interesting finding is the substantially differ- 
ent characterization of external shocks and entropy 
profiles in the grid and SPH methods, a feature that 
has a number of important consequences in both ther- 
mal and non thermal issues. The different dynamics 
felt by accreted clumps (Sec l5.4l) show that the pre- 
diction of mixing and gas matter deposition rates in 
cluster cosmological simulations is still an open prob- 
lem. Given the rather simple setup employed in these 
simulations (no radiative processes, no heating mech- 
anism other than shocks, no CR feedback, small tur- 
bulent motions due to lack of resolution and artificial 
viscosity in SPH) shocks dynamics has to be regarded 
as the leading player in setting the entropy profiles in 
clusters. These results conclusively suggest that the dif- 
ferences in shocks morphologies and shock dynamics 
across the clusters evolution leave major imprints also 
in substructures distributions and entropy distributions 
in the ICM, which is a rather new evidence provided 
by this work. 

Tightly connected to this, is the high degree of 
non-linearity which is present in all CR acceleration 
recipes. However, to date it is not clear whether these 
non-linearities would amplify any of the above differ- 
ences at shocks and potentially lead to a different CR 
pressure distribution in galaxy clusters, or whether the 
average CR pressure support results from a combina- 
tion of an average shock acceleration efficiency at the 
strongest shocks and successive CR transport. 

Based on the results of this project, we notice that 
performing of "high-precision" cosmology (e.g. relat- 
ing cosmological observables and theoretical models, 
based on scaling relations affected by less than ^ 1 
percent scatter in simulations) may still be a challenge 
for many applications, since some very important mea- 
surements related to the volume filling properties of 
galaxy clusters simulated with some of the best avail- 
able numerical codes appear to be still affected by un- 



certainties of ^ 10 percent (or more). This is found 
even in the case of the very similar physical setup an- 
alyzed here (only gravity and hydrodynamical forces) 
and the reason for this appear to be mostly of numerical 
nature, meaning that some important details concern- 
ing the production and transport of entropy in these 
simulations can be very different from code to code. 
This results are based on rather low or moderate resolu- 
tion simulations presented in this paper, and it is likely 
that going to much higher resolution levels off most of 
the above differences; however in this paper we have 
shown that not all significant differences are related to 
resolution only (e.g. differences in accretion shocks in 
SPH or grid methods). 

It is unclear if the application of more physical in- 
gredients which are not accounted in this work (e.g. 
magnetic fields, thermal conductions, feedback of rel- 
ativistic particles) may be able to soften any of the re- 
ported above reported differences. 

The suggestion of this work is that, together with 
the design of more sophisticated physical recipes to 
model the thermal and non thermal components of the 
real Universe, our theoretical understanding of cos- 
mic structures would also greatly benefit from other 
detailed comparative studies of different numerical 
recipes, since the convergence of simulated estimates 
of a sizable fraction of the cosmic volume is presently 
yet to be reached. 



8 ACKNOWLEDGEMENTS 

F.V. acknowledges support from the grant FOR 1254 
from Deutschen Forschungsgemeinschaft. K.D. ac- 
knowledges the support by the DFG Priority Pro- 
gramme 1 177 and additional support by the DFG Clus- 
ter of Excellence "Origin and Structure of the Uni- 
verse". The work of D. R. and H. K. was supported 
by the National Research Foundation of Korea through 
grant 2007-0093860. F.V. and G.B. acknowledge par- 
tial support through grant PRIN INAF 20082009/2010. 
We acknowledge the usage of computational resources 
under the CINECA-INAF 2008-2010 agreement and 
the 2009 Key Project "Turbulence, shocks and cosmic 
rays electrons in massive galaxy clusters at high resolu- 
tion". We acknowledge R.Brunino, M. Hoeft, T. Jones, 
B. O'Shea and V. Springel of fruitful scientific discus- 
sions. The computation of cosmological distances in 
this paper are done thanks to the Cosmology Calcula- 
tor by E. Wright (Wright 2006). 



REFERENCES 

Ackermann M. et al. , 2010, ApJL, 717, L71 
Agertz O., et al. , 2007, MNRAS, 380, 963 
Aharonian F., et al. , 2009, ApJ, 691, 175 
Aleksic, J., et al. 2010, ApJ, 710, 634 
Battaglia, N., Pfrommer, C., Sievers, J. L., Bond, J. R., 
& EnBlin, T. A. 2009, MNRAS, 393, 1073 



A Comparison of Cosmological Codes 29 



Berezinsky V. S., Blasi P., Ptuskin V. S., 1997, ApJ, 
487, 529 

Berger M. J., Colella P., 1989, JCoPh, 82, 64 
Berrington R. C, Dermer C. D., 2003, ApJ, 594, 709 
Blandford R. D., Ostriker J. P, 1978, ApJ, 221, L29 
Blasi P, Colafrancesco S., 1999, APh, 12, 169 
Blasi P, Gabici S., Brunetti G., 2007, IJMPA, 22, 681 
Borgani, S., Fabjan, D., Tomatore, L., Schindler, S., 
Dolag, K., & Diaferio, A. 2008, Space Science Re- 
views, 134, 379 

Borgani, S., & Kravtsov, A. 2009, arXiv:0906.4370 
Branetti G., Venturi T., Dallacasa D., Cassano R., 
Dolag K., Giacintucci S., Setti G., 2007, ApJ, 670, 
L5 

Brunetti G., et al. , 2008, Natur, 455, 944 
Brunetti, G., & Lazarian, A. 201 1, MNRAS, 410, 127 
Bryan G. L., Norman M. L., 1998, ApJ, 495, 80 
Bryan G. L., Norman M. L., 1997, ASPC, 123, 363 
Bullock, J. S., Kolatt, T. S., Sigad, Y., Somerville, 
R. S., Kravtsov, A. V., Klypin, A. A., Primack, J. R., 
& Dekel, A. 2001, MNRAS, 321, 559 
Bums, J. O., Skillman, S. W., & O'Shea, B. W. 2010, 
ApJ, 721, 1105 

Churazov, E., Forman, W., Vikhlinin, A., Tremaine, 
S., Gerhard, O., & Jones, C. 2008, MNRAS, 388, 
1062 

Colella, P, & Glaz, H. M. 1985, Journal of Computa- 
tional Physics, 59, 264 

Davis, M., Efstathiou, G., Frenk, C. S., & White, 
S. D. M. 1985, ApJ, 292, 371 
Dolag K., EnBlin T. A., 2000, A&A, 362, 151 
Dolag K., Vazza F, Brunetti G., Tormen G., 2005, 
MNRAS, 364, 753 

Dolag, K., Borgani, S., Schindler, S., Diaferio, A., & 
Bykov, A. M. 2008, Space Science Reviews, 134, 229 
Donnert, J., Dolag, K., Brunetti, G., Cassano, R., & 
Bonafede, A. 2010, MNRAS, 401, 47 
Drury L. O., Voelk J. H., 1981, ApJ, 248, 344 
Efstathiou, G., & Eastwood, J. W. 1981, MNRAS, 
194, 503 

Efstathiou, G. 1985, MNRAS, 213, 29P 

Eke, V. R., Cole, S., Frenk, C. S., & Patrick Henry, J. 

1998, MNRAS, 298, 1145 

Ensslin T. A., Biermann P. L., Klein U., Kohle S., 
1998, A&A, 332, 395 

EnBlin T. A., Simon P., Biermann P. L., Klein U., 
Kohle S., Kronberg P P, Mack K.-H., 2001, ApJ, 549, 
39 

Ensslin T. A., Pfrommer C, Miniati F., Subramanian 
K., 2011, A&A, 527,99 

Ettori S., Dolag K., Borgani S., Murante G., 2006, 
MNRAS, 365, 1021 

Frenk, C. S., et al. 1999, ApJ , 525, 554 
Fryxell, B., et al. 2000, ApJs, 131, 273 
Godunov, S. K. 1959, Mat. Sbomik, 47, 271 
Haardt F, Madau P, 1996, ApJ, 461, 20 
Harten, A., & Hyman, J. M. 1983, Journal of Compu- 
tational Physics, 50, 235 

Heitmann, K., et al. 2008, Computational Science and 
Discovery, 1, 015003 

HeB, S., & Springel, V. 2010, MNRAS, 406, 2289 
Hockney, R. W., & Eastwood, J. W. 1981, Computer 



Simulation Using Particles, New York: McGraw-Hill, 
1981, 

Hoeft, M., Bruggen, M., Yepes, G., Gottlober, S., & 

Schwope, A. 2008, MNRAS, 391, 1511 

Kang H., Jones T. W., 2002, JKAS, 35, 159 

Kang H., Ryu D., Cen R., Ostriker J. P, 2007, ApJ, 

669, 729 

Kang, H., & Ryu, D. 2010, ApJ, 721, 886 

Kravtsov A. V., Nagai D., Vikhlinin A. A., 2005, ApJ, 

625, 588 

Landau L. D., Lifshitz E. M., 1966, hydr.book, 
Leveque, R. J. 1990, VKl, Computational Fluid Dy- 
namics, Volume 1 142 p (SEE N90-27989 22-34), 1, 
Loeb A., Waxman E., 2000, Natur, 405, 156 
Merlin, E., Buonomo, U., Grassi, T., Piovan, L., & 
Chiosi, C. 2010, A&A, 513, A36 
Miniati F, Ryu D., Kang H., Jones T. W., 2001, ApJ, 
559, 59 

Miniati F, Ryu D., Kang H., Jones T. W., Cen R., Os- 
triker J. P, 2000, ApJ, 542, 608 
Miniati F, Jones T. W., Kang H., Ryu D., 2001, ApJ, 
562, 233 

Miniati F, 2003, MNRAS, 342, 1009 
Mitchell, N. L., McCarthy, I. G., Bower, R. G., The- 
uns, T., & Grain, R. A. 2009, MNRAS, 395, 180 
Molnar, S. M., Ream, N., Haiman, Z., Bryan, G., 
Evrard, A. E., & Lake, G. 2009, ApJ, 696, 1640 
Nagai, D., & Lau, E. T. 2011, Apjl, 731, LIO 
Norman M. L., Bryan G. L., Harkness R., Bordner J., 
Reynolds D., O'Shea B., Wagner R., 2007, arXiv, 705, 
arXiv:0705.1556 

Norman, M. L. 2010, arXiv: 1005.1 100 

O'Shea B. W., Bryan G., Bordner J., Norman M. L., 

Abel T., Harkness R., Kritsuk A., 2004, astro, 

arXiv:astro-ph/0403044 

O'Shea B. W, Nagamine K., Springel V, Hemquist 

L., Norman M. L., 2005, ApJS, 160, 1 

Peebles, P J. E. 1978, A&A, 68, 345 

Petrosian, V, & Bykov, A. M. 2008, SSR, 134, 207 

Pfrommer C, EnBlin T. A., 2004, JKAS, 37, 455 

Pfrommer C, Springel V, EnBlin T. A., Jubelgas M., 

2006, MNRAS, 367, 113 

Pfrommer C, EnBlin T. A., Springel V, Jubelgas M., 
Dolag K., 2007, MNRAS, 378, 385 
Pfrommer C, EnBlin T. A., Springel V, 2008, MN- 
RAS, 385, 1211 

Pfrommer C, 2008, MNRAS, 385, 1242 
Pfrommer C, Jones T. W., 201 1, ApJ, 730, 22 
Pinzke A., Pfrommer C, 2010, MNRAS, 409, 449 
Pinzke A., Pfrommer C, Bergstrom L., 201 1, submit- 
ted, arXiv: 1105.3240 

Price, D. J. 2008, Journal of Computational Physics, 
227, 10040 

Powell, K. G., Roe, P. L., Linde, T. J., Gombosi, T. I., 

& de Zeeuw, D. L. 1999, Journal of Computational 

Physics, 154, 284 

Reimer O., 2004, JKAS, 37, 307 

Robertson, B. E., Kravtsov, A. V, Gnedin, N. Y, 

Abel, T., & Rudd, D. H. 2010, MNRAS, 401, 2463 

Roncarelli, M., Ettori, S., Dolag, K., Moscardini, L., 

Borgani, S., & Murante, G. 2006, MNRAS, 373, 1339 

Ryu, D., Ostriker, J. P, Kang, H., & Cen, R. 1993, 



30 F.Vazzay K.Dolag, D.Ryu, G.Brunetti, C.Gheller, H.Kang, C.Pfrommer 



ApJ, 414, 1 

Ryu D., Kang H., Hallman E., Jones T. W., 2003, ApJ, 
593, 599 

Sarazin C. L., 1999, ApJ, 520, 529 
Sheth R. K., Tormen G., 1999, MNRAS, 308, 119 
Simionescu, A., et al. 2011, Science, 331, 1576 
Skillman, S. W., O'Shea, B. W., Hallman, E. J., Bums, 
J. O., & Norman, M. L. 2008, ApJ, 689, 1063 
Skillman, S. W., Hallman, E. J., O'Shea, B. W., Bums, 
J. O., Smith, B. D., & Turk, M. J. 2011, Apj, 735, 96 
Spergel D. N., et al. , 2007, ApJS, 170, 377 
Springel, v., Yoshida, N., & White, S. D. M. 2001, 
NA, 6, 79 

Springel V., 2005, MNRAS, 364, 1105 

Springel V., 2010, MNRAS, 401, 791 

Springel, V. 2010, ARAA, 48, 391 

Steinmetz, M., & White, S. D. M. 1997, MNRAS, 

288, 545 

Tasker, E. J., Bmnino, R., Mitchell, N. L., Michielsen, 
D., Hopton, S., Pearce, F. R., Bryan, G. L., & Theuns, 
T. 2008, MNRAS, 390, 1267 

Urban, O., Wemer, N., Simionescu, A., Allen, S. W., 
Bohringer, H. 2011, MNRAS, 546 
Valdamini, R., Ghizzardi, S., & Bonometto, S. 1999, 
NewA, 4,71 

Valdamini, R. 2002, ApJ, 567, 741 

Vazza, R, Branetti, G., & Gheller, C. 2009, MNRAS, 

395, 1333 

Vazza, R, Bmnetti, G., Kritsuk, A., Wagner, R., 
Gheller, C., & Norman, M. 2009, A&A, 504, 33 
Vazza R, 201 1, MNRAS, 410, 461 
Vazza R, Bmnetti G., Gheller C., Bmnino R., 2010, 
NewA, 15, 695 

Vazza R, Gheller C, Bmnetti G., 2010, A&A, 513, 
A32 

Vazza, R, Roncarelli, M., Ettori, S., & Dolag, K. 
2011, MNRAS, 615 

Voit, G. M., Kay, S. T., & Bryan, G. L. 2005, MN- 
RAS, 364, 909 

Volk H. J., Aharonian F. A., Breitschwerdt D., 1996, 
SSRv, 75, 279 

Wadsley, J. W, Veeravalh, G., & Couchman, H. M. R 
2008, MNRAS, 387, 427 

Warren, M. S., Abazajian, K., Holz, D. E., & Teodoro, 
L. 2006, ApJ, 646, 881 

Woodward R, Colella R, 1984, JCoPh, 54, 1 15 



